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The  general  goal  of  this  research  was  to  model  the  translation  of  motor  cortical  commands  for 
movements  in  space  to  motoneuronal  output  through  an  intercalated  system  of  intercalated 
intemeurons. 

The  specific  objectives  were: 

1.  To  develop  a  physiologically  relevant,  neural  network  to  model  time- varying  neuronal 
population  operations  in  the  motor  cortex,  dealing  with  movements  in  space; 

2.  To  develop  a  model  of  interactions  between  the  cortical  and  the  spinal  networks  dealing 
with  generating  time- varying  motoneuronal  outputs  for  movements  in  space. 

The  novelty  of  our  approach  consisted  in: 

(a)  The  realistic  nature  of  the  elements  in  our  networks; 

(b)  The  massive  and  asymmetric  intercormectivity  among  network  elements; 

(c)  The  physiologically  relevant  design  of  the  networks,  including  the  communication  by 
spike  trains  among  network  elements  and  rules  of  cormectivity  based  on  experimental 
findings; 

(d)  The  dynamical  behavior  of  the  networks;  and 

(e)  The  time- varying  performance  of  the  networks. 

Results  and  Discussion 

1.  Development  of  a  physiologically  relevant,  neural  network  to  model  time-varying 
neuronal  population  operations  in  the  motor  cortex,  dealing  with  movements  in  space 
(publications  1-6;  numbers  in  parentheses  refer  to  the  list  of  papers  in  which  the  work  was 
published,  whereas  named  quotations  refer  to  the  Literature  Cited). 

We  published  7  papers  on  this  work.  We  summarize  below  the  main  findings  with  respect  to  (A) 
a  network  using  simple  input-output  neurons,  (B)  a  network  using  biophysically  modeled, 
spiking  neurons,  and  (C)  a  model  of  interacting  neural  networks. 
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We  developed  physiologically  relevant,  neural  networks  to  model  time-varying  neuronal  population 
operations  in  the  motor  cortex  and  spinal  cord,  dealing  with  movements  in  space.  We  also  developed 
a  model  of  the  interactions  between  these  two  networks  dealing  with  generating  time-varying 
motoneuronal  outputs  for  movements  in  space.  The  novelty  of  our  approach  consisted  in  (a)  the  realistic 
nature  of  the  elements  in  our  networks,  (b)  the  massive  and  asymmetric  intercoimectivity  among 
network  elements,  (c)  the  physiologically  relevant  design  of  the  networks,  including  the  communication 
by  spike  trains  among  network  elements  and  rules  of  connectivity  based  on  experimental  findings,  (d) 
the  dynamical  behavior  of  the  networks,  and  (e)  the  time-varying  performance  of  the  networks.  Finally, 
we  were  able  to  reliably  decode  and  transform  the  neuronal  ensemble  activity  recorded  in  behaving 
animals  for  controlling  an  simulated  arm.  This  demonstration  suggests  that  the  use  of  biologically 
inspired  neural  networks  to  transform  raw  cortical  signals  into  the  motor  output  of  a  multijoint  artificial 
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I.A.  Motor  cortical  network  using  simple  input-output  neurons 


Abstract 

As  a  dynamical  model  for  motor  cortical  activity  during  hand  movement  we  constructed  an 
artificial  neural  network  that  consists  of  extensively  interconnected  neuron-like  units  and 
performs  the  neuronal  population  vector  operations.  Local  geometrical  parameters  of  a  desired 
curve  are  introduced  into  the  network  as  an  external  input.  The  output  of  the  model  is  a  time- 
dependent  direction  and  length  of  the  neuronal  population  vector  which  is  calculated  as  a  sum  of 
the  activity  of  directionally  tuned  neurons  in  the  ensemble.  The  main  feature  of  the  model  is  that 
dynamical  behavior  of  the  neuronal  population  vector  is  the  result  of  connections  between 
directionally  tuned  neurons  rather  than  being  imposed  externally.  The  dynamics  is  governed  by  a 
system  of  coupled  non-linear  differential  equations.  Connections  between  neurons  are  assigned 
in  the  simplest  and  most  common  way  so  as  to  fulfill  basic  requirements  stemming  fi:om 
experimental  findings  concerning  the  directional  tuning  of  individual  neurons  and  the 
stabilization  of  the  neuronal  population  vector,  as  well  as  fi-om  previous  theoretical  studies.  The 
dynamical  behavior  of  the  model  reveals  a  close  similarity  with  the  experimentally  observed 
dynamics  of  the  neuronal  population  vector.  Specifically,  in  the  firamework  of  the  model  it  is 
possible  to  describe  a  geometrical  curve  in  terms  of  the  time  series  of  the  population  vector.  A 
correlation  between  the  dynamical  behavior  of  the  direction  and  the  length  of  the  population 
vector  entails  a  dependence  of  the  "neural  velocity"  on  the  curvature  of  the  tracing  trajectory  that 
corresponds  well  to  the  experimentally  measured  covariation  between  tangential  velocity  and 
curvature  in  drawing  tasks. 

Introduction 

Although  individual  neurons  in  the  motor  cortex  are  broadly  tuned  with  respect  to  the  direction 
of  movement  (Georgopoulos  et  al.  1982)  and  a  particular  direction  involves  the  activation  of  a 
large  ensemble  of  cells,  the  high  accuracy  of  a  movement  performance  can  be  explained  by 
means  of  "population  coding"  hypothesis  (Georgopoulos  et  al.  1983;  Georgopoulos  et  al.  1986; 
Georgopoulos  et  al.  1988).  In  this  hypothesis  the  contribution  of  each  i'*  directionally  timed 
neuron  is  represented  as  a  vector  pointing  in  its  preferred  direction  C,-  with  the  length  given  by 
the  change  in  cell  frequency  of  discharge  associated  with  a  particular  movement  direction.  For  a 
given  movement  M,  these  contributions  add  vectorially  to  yield  the  "neuronal  population  vector" 
P,  which  is  a  measure  of  the  combined  directional  tendency  of  the  whole  neuronal  ensemble: 


P(M;t)  =  y;v,(M;t)C, 


(1) 


where  F^(M;r)  is  the  activity  of  the  f  neuron  at  time  bin  t  (usually  10-ms  or  20-ms  bin).  Now  it  is 
well-established  that  the  neuronal  population  vector  predicts  accurately  the  direction  of 
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movement  under  a  variety  of  conditions  (Georgopoulos  1990). 

Remarkably,  the  population  vector  can  be  used  as  a  probe  to  monitor  the  directional  tendency 
of  the  neuronal  ensemble,  as  it  changes  in  time.  For  example,  in  a  task  that  required  a  movement 
to  be  made  at  a  given  angle  from  a  stimulus  direction,  the  population  vector  rotated  from  the 
stimulus  direction  to  the  direction  of  the  movement  (Georgopoulos  et  al.  1989;  Lurito  et  al. 

1991).  Moreover  the  "neural  trajectories"  reconstructed  from  the  time  series  of  neuronal 
population  vectors  attached  to  each  other  tip-to-tail  correspond  well  to  the  actual  movement 
trajectories  produced  in  reaching  (Georgopoulos  et  al.  1988)  as  well  as  in  drawing  tasks 
(Schwartz  1993, 1994).  All  of  these  findings  indicate  that  the  ongoing  direction  and  length  of  the 
neuronal  population  vector  could  predict  the  tangential  velocity  vector  of  the  movement  that  is 
made  approximately  120  ms  later  (Schwartz  1993).  In  its  turn,  a  stable  covariation  between  the 
tangential  velocity  of  the  hand  movement  and  the  curvature  of  a  drawn  trajectory  has  been  well 
documented  (Viviani  and  Terzuolo  1982;  Lacquaniti  et  al.  1983;  Soechting  and  Terzuolo  1986; 
Massey  et  al.  1992).  Qualitatively,  this  correlation  reflects  the  fact  that  the  movement  tends  to 
slow  down  when  the  trajectory  is  more  curved.  It  was  hypothesized  (Massey  et  al.  1992)  that  the 
underlying  mechanisms  of  the  covariation  between  geometrical  and  kinematic  parameters  of  the 
movement  could  be  neural  constraints.  Namely,  if  indeed  the  neuronal  population  vector  predicts 
the  tangential  velocity  vector  of  the  real  movement,  then  the  reason  of  the  tendency  to  slow  down 
at  the  curved  sections  of  a  trajectory  could  be  explained  by  the  additional  time,  compared  with 
straight-line  movement,  needed  for  the  rotation  of  the  neuronal  population  vector.  However,  the 
neural  basis  which  would  explain  qualitative  and  quantitative  aspects  of  the  problem  remains  to 
be  explored. 

The  distributed  coding  of  the  direction  of  movement  as  well  as  stationary  properties  of  the 
ensemble  of  motor  cortical  units  have  been  well  reproduced  in  the  framework  of  simple  three¬ 
layered  model  trained  through  a  supervised  learning  algorithm  (Lukashin  1990),  by  means  of  the 
model  of  arrays  of  adjustable  pattern  generators  (Berthier  et  al.  1991)  and  within  a  more 
biological  model  that  computes  visuomotor  transformations  for  arm  reaching  movements 
(Bumod  et  al.  1992).  The  purpose  of  this  paper  is  to  propose  a  simple  dynamical  neural  network 
model  which  consists  of  extensively  interconnected  neuron-like  units  and  performs  the  neirronal 
population  vector  operations  as  a  response  to  some  input  information.  Specifically,  we  consider  a 
model  of  neural  network  which  receives  the  input  in  terms  of  geometrical  features  of  a  trajectory 
and  produces  the  dynamical  behavior  of  the  neuronal  population  vector  as  the  output.  The 
dynamical  properties  of  the  neuronal  population  vector  as  well  as  a  relation  between  the 
curvature  of  a  trajectory  and  tangential  velocity  obtained  in  the  framework  of  the  model  reveal 
close  correspondence  with  the  experimental  findings  mentioned  above. 


Model 


Basic  requirements.  Generally,  the  present  model  belongs  to  a  well-known  class  of  neural 
networks  which  are  able  to  generate  sequential  patterns  of  neural  activity  (Kleinfeld  1986; 
Sompolinsky  and  Kanter  1986;  Jordan  1986;  Dehaene  et  al.  1987;  Kleinfeld  and  Sompolinsky 
1988;  Guyon  et  al.  1988;  Schoner  and  Kelso  1988;  Massone  and  Bizzi  1989;  Williams  and 
Zipser  1989;  Pearlmutter  1989;  Amirikian  and  Lukashin  1992).  Therefore  below  we  will 
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concentrate  only  on  those  features  of  the  model  that  are  essential  for  the  problem  in  question. 

Let  and  Ky  be  the  x  and  y  components  of  a  unit  length  tangential  vector  K  taken  at  a  point 
K  on  drawn  trajectory.  Let  and  My  be  the  same  values  for  a  subsequent  point  M.  Receiving 

components  of  the  vectors  K  and  M  as  external  input,  the  desired  network  has  to  produce 
components  of  the  neuronal  population  vector  P  as  a  sum  of  the  responses  of  the  ensemble 
directionally  tuned  neurons  [see  (1)].  The  basic  qualitative  requirements  for  the  dynamical 
behavior  of  the  neuronal  population  vector  follow  from  experimental  findings  concerning  the 
dynamical  behavior  of  the  neuronal  population  vector  during  a  transformation  task 
(Georgopoulos  et  al.  1989;  Lurito  et  al.  1991).  The  observed  dynamical  behavior  of  the  neuronal 
population  vector  in  the  case  with  externally  assigned  initial  (stimulus)  and  final  (movement) 
directions  consists  of  the  rotation  of  the  population  vector  from  an  initial  to  a  final  direction, 
followed  by  the  stabilization  of  the  direction  and  the  length  of  the  population  vector.  Hence  we 
suppose  that  at  the  motor  cortical  level  the  tracing  out  the  section  of  a  trajectory  between  point  K 
with  tangential  vector  K  and  point  M  with  tangential  vector  M  corresponds  to  the  rotation  of  the 
neuronal  population  vector  P  from  the  direction  K  to  the  direction  M  with  subsequent 
stabilization  in  the  direction  M.  The  stabilization  of  the  direction  of  the  population  vector  triggers 
the  change  of  the  external  signals:  now  the  direction  M  is  considered  by  the  network  as  the  initial 
direction  and  a  new  final  desired  direction  (for  example,  tangential  vector  N  of  the  next  point  N 
on  the  trajectory)  is  given  as  a  new  external  input.  This  assumption  is  conceptually  in  accord 
with  the  mechanisms  of  continuous  updating  of  the  current  position  commands  to  arm  system 
muscles,  introduced  by  Bullock  and  Grossberg  (1988)  in  a  model  of  the  circuit  that  automatically 
generates  arm  movement  trajectories. 

Architecture.  The  input  (upper)  layer  of  the  network  consists  of  four  neurons  the  activities 
of  which  code  components  of  vectors  K  and  M.  Only  the  activities  of  these  neurons  are  assigned 
externally.  The  output  (lowest)  layer  consists  of  two  neurons  with  linear  activation  fimctions. 

The  activities  of  these  neurons  are  supposed  to  code  components  of  the  neuronal  population 
vector  P.  There  are  two  hidden  layers  in  the  network:  a  layer  of  modulators  whose  activities  are 
designated  by  Q  (only  one  of  them  is  shown)  and  a  layer  of  directionally  tuned  neurons  whose 
activities  are  designated  by  V^,  where  index  i  enumerates  neurons  in  the  layer.  The  activities  of  all 
neurons  in  both  hidden  layers  are  determined  by  a  nonlinear  activation  function  which  is  chosen 
as  hyperbolic  tangent  function. 

There  are  connections  v^,  v,y  between  neurons  Ky  of  the  input  layer,  which  code  the  initial 
direction  K,  and  all  neurons  Vi  of  directionally  tuned  hidden  layer.  Correspondingly,  there  are 
connections  Wy,  between  all  Fj  neurons  and  the  two  neurons  Py  in  the  output  layer.  The 
essential  feature  of  the  model  is  the  existence  of  intralayer  connections  between  neurons  F,  in  the 
directionally  tuned  layer.  These  connections  are  of  two  types:  modulated,  Vy,  and  non-modulated, 
Wy  (see  below).  The  modulator  neuron  Q  receives  inputs  from  neurons  M,,  My,  which  code  the 
required  final  direction  M,  and,  by  means  of  feedback  connections,  from  the  two  neurons  P^^,  Py 
of  the  output  layer,  which  code  the  current  direction  of  the  neuronal  population  vector.  For  the 
sake  of  simplicity  the  absolute  values  of  all  these  coimections  are  supposed  to  have  the  same 
strength  q'"'’.  The  modulator  neuron  Q  sends  signals  to  directionally  tuned  neurons  F^-  in  a  way 
that  the  efficacy  of  modulated  connections  Vy  between  neurons  Fj-  and  the  efficacy  of  connections 
Vfe,  Vjy  between  neurons  F,  and  neurons  Ky  depend  on  the  activity  of  modulator  Q  and  on  the 
connection  strengths  between  neuron  Q  and  neurons  F).  Again  the  latter  connections  are 
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supposed  to  have  the  same  strength 

Dynamical  equations.  In  the  framework  of  standard  theory  of  neural  networks  (see,  for 
instance,  Amari  1972;  Sejnowski  1976;  Grossberg  and  Cohen  1983;  Hopfield  1984;  Hopfield 
and  Tank  1986;  Atiya  and  Baldi  1989)  the  dynamical  behavior  of  the  model  described  above  is 
governed  by  the  following  system  of  coupled  nonlinear  differential  equations: 


^--Ui(t)  +  ^WijVj(t)  + 

at  j 


2;VijVj(t)  +  VixKx  +  ViyKy 


V  J 


q°"‘Q(t)  (2) 


F,(if)  =  tanh(ui(t)) 


(3) 


Px(if)  =  SwxiVi(t);  Py(t)  =  XwyiVi(t)  (4) 

i  i 

Q(t)  =  tanhlq'"”  (  [Mx  -  Px  (t)|  +  |My  -  Py  (t)| )]  (5) 


Argximent  t  is  shown  for  variables  which  depend  on  time.  Equation  (2)  is  a  resistance-capacitance 
equation  for  the  z'*  neuron  of  the  directionally  tuned  layer  (variable  for  example,  might 
represent  the  soma  membrane  potential  of  the  neuron  averaged  over  a  reasonable  time  interval). 
The  time  constant  <5  is  a  characteristic  time  for  the  individual  neuron.  To  simplify  the  equations, 
characteristic  times  for  neurons  in  other  layers  are  supposed  to  be  much  shorter  then  the  6  value. 
We  checked  that  this  approximation  does  not  change  die  main  results  presented  below.  Equation 
(3)  shows  the  relation  between  internal  state  of  the  i*’'  neuron  and  its  output  activity  (this  value 
might  correspond  to  the  change  of  discharge  rate  of  a  real  cell  averaged  over  a  reasonable  time 
interval).  In  accordance  with  the  model,  the  activity  of  the  modulator  Q  at  time  t  (5)  is  the 
function  of  the  components  of  the  constant  vector  M  and  the  components  of  vector  P  at  time  t. 
[Tilde  above  means  that  vector  P  is  normalized  to  unity  in  (4)].  The  components  of  vector  P  are 
linear  functions  of  activities  F((/),  as  can  be  seen  from  (4). 

Once  equations  (2)-(5)  are  written,  the  dynamical  behavior  of  the  model  completely  depends 
on  the  set  of  connection  strengths.  In  the  next  section  we  explain  the  motivation  for  our  choice  of 
these  parameters. 


6 


Connection  strengths 

Directional  tuning  of  neurons.  Suppose  that  there  are  no  intralayer  connections  between 
neurons  (or  Wy  =  0,  Vy  =  0)  and  no  feedback  connections  to  modulator  (Q(t) .  1).  Then 
equations  (2)  are  reduced  to  ordinary  differential  equations  having  straightforward  solution, 
which  for  r  o  d  and  q°'“  n  1  leads  to  the  relations: 


Vi(t  «T ;  «  {vuK,  +  ViyKy)=  q°“‘ki  cos(e k -ai) 


(6) 


where  e^^-  is  the  angle  characterizing  the  direction  of  the  input  vector  K,  and  parameters  and  a,- 
are  expressed  through  connection  strengths  v,y.  Since  (6)  gives  the  same  response  activity  to 
the  external  direction  K  as  has  been  observed  in  numerous  experiments  (Georgopoulos  et  al. 
1982;  Schwartz  et  al.  1988;  Schwartz  1992)  the  angle  a,  can  be  regarded  as  the  preferred 
direction  of  the  i“'  neuron.  The  important  property  of  the  real  ensemble  of  directionally  tuned 
neurons  is  the  uniform  distribution  of  preferred  directions  in  space  (Schwartz  et  al.  1988;  Bumod 
et  al.  1992).  The  same  type  of  distribution  was  obtained  in  the  framework  of  the  simplified 
model  (without  intralayer  connections)  by  means  of  adjusting  the  set  of  the  parameters 

through  a  supervised  learning  algorithm  (Lukashin  1990).  In  accordance  with  these  results,  we 
use  for  the  complete  model  (2)-(5)  the  following  expressions  for  connection  strengths  [compare 
with  (6)]: 


V£,-cosai;  Viy-sinai 


(7) 


where  angles  a,  are  randomly  and  uniformly  distributed  over  the  interval  [-8,6].  Since  the 
preferred  directions  are  given,  connection  strengths  w^,-,  (see  (4))  should  lead  to  the  standard 
definition  of  the  neuronal  population  vector  (1).  Namely, 


2  2 

>v«  =  — costti;  Wyi  =  — sinai 
N  N 


(8) 


where  (2/7V)  is  a  normalizing  coefficient,  N  is  the  number  of  neurons  in  the  directionally  tuned 
layer.  Equations  (8),(4)  are  completely  equivalent  to  (1). 

Stabilization  of  the  neuronal  population  vector.  Suppose  that  the  direction  of  the  neuronal 
population  vector  P  is  close  enough  to  the  desired  final  direction  Af,  or  Q(t)  .  0  [see  (4)].  Then 
only  the  first  and  second  terms  in  the  right-hand  side  of  (2)  are  retained.  In  accordance  with 
experimental  findings  in  this  case  (Georgopoulos  et  al.  1989;  Lurito  et  al.  1991),  the  direction  of 
the  neuronal  population  vector  should  be  stable.  During  this  steady-state  period  activities  of 
directionally  neurons  cease  to  change,  so  that  du/dt  =  0  for  all  i,  which  together  with  (2)  and  (3) 
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leads  to  the  condition 


F,YM)  =  tanh 


SwijVjCM) 


V  j 


(9) 


that  should  be  valid  for  all  neurons  and,  generally,  for  all  directions  M.  If  activities  F)-  at  the 
stable  state  are  known  for  a  given  number  of  neurons  and  for  an  appropriate  number  of  directions 
M,  equation  (9)  can  be  used  for  determining  the  parameters  w^.  The  general  features  of  the  sets  of 
connection  strengths  Wy  that  would  ensure  the  stability  was  extensively  studied  (Georgopoulos  et 
al.  1993).  It  was  shown  that  imder  some  restrictions  on  its  possible  values,  the  Wy  parameters 
ensuring  condition  (9)  depend  on  preferred  directions  of  the  pair  of  neurons  involved  in  the 
connection.  Namely,  the  mean  synaptic  strength  was  negatively  correlated  with  the  angle 
between  the  preferred  directions  of  the  two  neurons  throughout  the  range  of  connections,  from 
positive  (excitation)  to  negative  (inhibition).  These  results  of  modeling  were  confirmed  by  the 
findings  of  neurophysiological  studies  [for  details  see  (Georgopoulos  et  al.  1993)].  The  simplest 
way  to  include  the  dependence  on  preferred  directions  is  to  approximate  the  Wy  parameters  by  a 
cosine  relation: 


Wj,--  — cos(ai-aj) 


(10) 


However  we  checked  that  sets  of  the  parameters  Wy  obtained  in  (Georgopoulos  et  al.  1993)  lead 
to  the  same  results  as  the  set  (10). 

Rotation  of  the  neuronal  population  vector.  The  dynamical  behavior  of  the  system 
producing  the  rotation  of  the  neuronal  population  vector  can  be  ensured  by  connections  that 
induce  the  phase  shift  violating  steady-state  conditions.  Again  the  simplest  way  to  include  this 
phase  shift  is  to  assign  modulated  intralayer  connection  strengths  Vy  by  the  relations: 


v//  =  — sin(ai-aj),  (i_j) 


(11) 


To  allow  for  the  possibility  of  changing  the  length  of  the  population  vector  we  introduce  a  self¬ 
inhibition  for  activities  of  neurons  in  the  directionally  tuned  layer  (non-zero  negative  values  of 
the  diagonal  elements  v„).  The  self-inhibition  is  under  control  by  the  feedback  loop  in  the  same 
way  as  it  is  realized  for  the  modulator  Q . 


Results  of  simulations  and  discussion 
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Below  we  give  the  results  of  calculations  of  the  dynamical  behavior  of  the  model  described  by 
(2)-(5)  with  connection  strengths  given  by  (7),(8),(10),(1 1).  Preferred  directions  a,-  are  randomly 
and  uniformly  distributed  over  the  interval  [-6,6].  The  number  of  neurons  in  directionally  tuned 
layer  in  routine  calculations  is  A/’  =  1 00. 

Adjustable  parameters.  Only  three  adjustable  parameters  remain  to  be  established:  the 
characteristic  time  6  (2),  the  input,  9'"'’,  and  the  output,  q°“',  connections  of  the  modulator  [see 
(2),(5)].  The  parameter  6  determines  the  time  scale  of  the  dynamics.  A  reasonable  value  is  d .  5 
ms.  We  checked  that  the  results  do  not  depend  on  a  particular  value  of  the  q'"’’  parameter  if  it  is 
large  enough.  (All  the  results  presented  below  have  been  obtained  for  q'”^  =  50.)  In  contrast,  the 
model  reveals  a  strong  sensitivity  to  the  value  of  the  q°'“  parameter.  The  chosen  set  of  connection 
strengths  provides  dynamics  such  as  a  movement  towards  an  attractor.  It  can  be  seen  that  the 
population  vector  rotates  from  the  direction  of  vector  K  towards  the  direction  of  vector  M  and 
then  stabilizes  in  the  latter  direction.  The  shape  of  the  curves  as  well  as  the  saturation  level 
essentially  depends  on  the  value  of  the  q°‘‘^  parameter.  It  makes  sense  to  adjust  this  parameter  so 
as  to  obtain  the  angular  velocity  of  rotation  at  initial  moments  of  time  close  to  the  experimentally 
measured  value  which  is  of  the  order  of  500  deg/s  (Georgopoulos  et  al.  1989;  Lurito  et  al.  1991). 
From  this  point  of  view  the  curve  corresponding  to  q'^'“  =  0.05  gives  the  best  fit.  Indeed,  in  this 
case  the  angular  velocity  is  around  2.5  deg  per  6  unit  or  500  deg/s,  if  d  .  5  ms.  All  the  results 
shown  below  were  obtained  for  q°'“  =  0.05. 

Feedback  connections.  In  the  beginning  of  the  dynamics,  when  the  current  direction  of 
population  vector  is  far  enough  from  the  desired  final  direction,  the  population  vector  rotates 
with  constant  angular  velocity.  When  the  direction  of  population  vector  becomes  close  to  the 
assigned  final  direction  M,  the  activity  of  the  modulator  Q  tends  to  zero.  The  third  term  in  the 
right-hand  side  of  the  basic  equation  (2)  ceases  to  affect  the  dynamics,  and  the  process  enters  into 
the  steady-state  phase.  The  stabilized  directions  of  the  population  vector  coincide  with  the 
directions  of  vectors  M.  In  the  light  of  experimental  data  concerning  the  dynamical  behavior  of 
population  vector  during  a  transformation  task  (Georgopoulos  et  al.  1989)  we  consider  the 
steady-state  phase  as  the  necessary  part  of  the  whole  dynamics. 

Population  coding  of  movement  trajectories.  The  aim  of  this  section  is  to  check  the 
ability  of  the  model  to  describe  trajectories  in  terms  of  time  series  of  the  neuronal  population 
vectors  attached  to  each  other  tip-to-tail  if  the  geometrical  features  of  a  trajectoiy  are  introduced 
into  the  network  as  external  input.  Here,  as  an  example,  we  consider  closed  elliptic  trajectories. 
We  checked  that  all  the  results  remain  valid  for  other  curves. 

Let  a  desired  trajectory  be  given.  Let  K,  M,  N,...  be  points  along  the  trajectory  separated  from 
each  other  by  equal  arc  lengths  (in  routine  calculations  we  used  100  points  per  each  trajectory). 
The  beginning  of  the  dynamics  is  considered  to  be  the  instant  at  which  the  components  of  first 
two  tangential  vectors  K  (initial)  and  M  (final)  are  introduced  into  the  network  as  activities  of 
corresponding  neurons.  Then  the  dynamics  described  by  (2)-(5),(7),(8),(10),(l  1)  begins.  When 
the  current  direction  of  population  vector  P  becomes  close  to  the  direction  M  the  activity  of  the 
modulator  Q  rapidly  tends  to  zero.  On  the  one  hand,  this  entails  the  stabilization  of  the 
population  vector  P.  On  the  other  hand,  the  low  value  of  the  activity  of  the  modulator  Q  triggers 
the  change  of  the  input  signals.  Now  vector  M  serves  as  the  initial  direction,  and  the  next 
tangential  vector  N  is  introduced  as  a  desired  final  direction,  and  so  on.  It  is  important  to  note 
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that  in  the  framework  of  this  procedure  the  relative  velocity  of  tracing  of  the  trajectory  is 
determined  by  the  network  itself  because  the  instant  of  appearance  of  the  next  input  signal 
depends  only  on  the  internal  state  of  the  network. 

Conclusion 

The  main  feature  of  the  proposed  model  is  that  the  neuronal  population  vector  operations 
produced  by  the  network  are  the  result  of  the  interactions  between  directionally  tuned  neurons 
rather  than  being  imposed  by  external  inputs.  Connections  between  neurons  were  assigned  in  the 
simplest  and  most  common  way  so  that  actually  only  one  numerical  parameter,  namely  has 
been  adjusted.  The  network  receives  inputs  as  a  sequence  of  tangential  vectors  of  a  curve  which 
is  to  be  traced  out,  but  the  "neural  velocity"  of  the  tracing  is  a  function  of  internal  parameters  of 
the  model.  The  output  of  the  network  is  the  current  direction  and  length  of  the  population  vector. 
The  correspondence  between  the  changes  of  these  parameters  during  tracing  a  trajectory  makes  it 
possible  to  describe  a  drawn  geometrical  curve  in  terms  of  the  time  series  of  the  population 
vector  in  full  agreement  with  experimental  findings  (Georgopoulos  et  al.  1988;  Schwartz  1993, 
1994).  Moreover,  the  dependence  of  the  "neural  velocity"  on  the  curvature  obtained  in  the 
framework  of  the  model  reveals  close  correspondence  with  the  similar  covariation  observed 
experimentally  for  hand  movement  in  drawing  tasks  (Viviani  and  Terzuolo  1982;  Lacquaniti  et 
al.  1983;  Soechting  and  Terzuolo  1986;  Massey  et  al.  1992). 

The  network  considered  in  this  paper  possesses  three  major  characteristics  that  are 
biologically  relevant  for  cortical  operations  in  general  and  two  properties  that  are  important  for 
the  motor  cortex  in  particular.  The  general  characteristics  include  the  massive  intercormectivity 
among  the  network  elements,  the  weak  strength  of  synaptic  coimections  and  the  correlated 
interactions:  indeed,  these  seems  to  be  general  features  of  cortical  processing  (Martin  1988).  The 
specific  properties  of  the  network  that  correspond  to  motor  cortical  operations  include  the 
directional  tuning  of  individual  elements  and  the  computation  of  the  neuronal  population  vector. 
Altogether  then,  this  network  combines  biologically  meaningful  aspects  which  are  embedded 
within  a  dynamical  behavior  that  evolves  in  time.  This  temporal  evolution  combines  an 
experimentally  observed  operation  (i.e.  the  rotation  of  the  neuronal  population  vector)  and 
applies  it  to  the  neural  representation  of  curvature  in  a  planned  movement  trajectory.  However, 
for  this  application  an  internal  feedback  mechanism  is  proposed  that  signals  the  completion  of  a 
unit  rotation;  this  hypothesized  feedback  mechanism  needs  to  be  identified  experimentally.  It  is 
this  interplay  between  theory  and  experiment  that  lies  in  the  heart  of  the  present  attempt  to  bridge 
artificial  and  real  neural  networks  within  the  domain  of  the  motor  cortex. 


l.B.  Motor  cortical  network  using  biophysically  modeled,  spiking  neurons 

Abstract 

A  biophysically  realistic  neural  network  is  presented  that  models  operations  in  an  ensemble  of 
directionally  tuned  neurons  in  the  motor  cortex.  The  model  reproduces  well  directional 
operations  previously  identified  experimentally,  including  the  prediction  of  the  direction  of  an 


10 


upcoming  movement  in  reaching  tasks  and  the  rotation  of  the  neuronal  population  vector  in  a 
directional  transformation  task. 

Introduction 

The  model  described  in  the  previous  section  consisted  of  formal,  input-output  neurons.  In  this 
work  we  describe  a  biologically  relevant  model  neural  network  that  models  the  motor  cortical 
activity  during  the  intended  movement  in  the  framework  of  the  second  hypothesis  above.  The 
key  idea  is  that  a  correlation  between  synaptic  connection  strengths  and  preferred  directions  of 
the  neurons  involved  in  the  connection  could  ensure  the  observed  stability  and  transformations  of 
the  neuronal  population  vector.  We  showed  previously  (Georgopoulos  et  al.  1993)  that  this 
holds  for  a  formal  neural  network;  in  this  paper  we  extend  this  work  to  a  realistically  modeled, 
spiking  neural  network. 

Model 

Although  simulations  of  the  dynamical  behavior  of  the  neuronal  population  vector  by  modeling 
neurons  as  simple  input-output  devices  were  carried  out  (Lukashin  1990;  Lukashin  and 
Georgopoulos  1993),  it  is  still  unclear  where  is  the  border  between  important  features  and 
incidental  details  in  the  modeling  of  motor  cortical  cells.  Hence,  more  realistic  models  are 
desirable.  In  the  present  paper  we  simulate  an  artificial  network  that  consists  of  48  heavily 
interconnected  directionally  tuned  neurons  modeled  at  a  reasonable  level  of  reality.  Real  motor 
cortical  neurons  produce  trains  of  action  potentials,  or  spikes,  and  it  is  the  timing  of  these  spikes 
that  carries  information  about  the  directional  tendency.  At  the  same  time,  for  the  problem  in 
question  it  is  obviously  premature  to  allow  for  details  concerning,  for  example,  spatial 
localization  of  channels  or  dendritic  branching  patterns. 

To  describe  a  neuron  and  interactions  between  cells  we  use  a  biophysical  model  of 
intermediate  complexity.  A  neuron  is  represented  as  a  one-compartment  cell  equipped  with  a 
standard  set  of  active  chaimels  with  kinetics  governed  by  the  Hodgkin-Huxley  type  equations 
(Hodgkin  and  Huxley  1952;  Koch  and  Segev  1989).  At  present,  each  cell  has  voltage-dependent 
sodium  and  potassium  channels,  and  a  calcium-dependent  potassium  channel  that  generates  the 
late-phase  afterhyperpolarization.  The  intracellular  level  of  free  calcium  is  the  sum  of  two 
calcium  currents:  the  first  corresponds  to  the  influx  through  voltage  gated  calcium  channel,  and 
the  second  corresponds  to  calcium  that  enters  through  the  NMDA  channel.  The  depolarizing 
current  resulting  from  the  influx  of  calcium  ions  has  been  ignored.  Specifically,  the  system  of 
equations  listed  in  (Ekeberg  et  al.  1991)  has  been  used  to  calculate  the  ionic  currents,  calcium 
level  and  corresponding  voltage-  and  time-  dependent  degrees  of  activation/inactivation  of  the 
channels.  The  relevant  parameters  are  presented  in  Table  1.  (The  28  parameters  were  used  in 
describing  exponential  expressions  for  the  rate  constants  of  activation/inactivation  of  the 
chaimels  are  exactly  the  same  as  those  listed  in  Table  1  of  (Ekeberg  et  al.  1991)  and  are  not 
shown). 

All  neurons  are  connected  with  each  other,  and  changes  in  synaptic  conductance  are  used 
to  model  conventional  excitatory  and  inhibitory  synaptic  interactions.  A  contribution  of  the 
synaptic  currents  into  the  total  current  entering  a  cell  is  calculated  as  a  liner  sum  over  all 
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synapses.  Ionic  current  resulting  from  activation  of  a  synapse  due  to  arriving  of  a  presynaptic 
action  potential  is  calculated  as  the  product  of  time-dependent  synaptic  conductance  and  the 
difference  between  the  reversal  potential  of  the  synapse  and  the  membrane  potential  of  the  cell. 
The  synaptic  conductance  decreases  exponentially  with  a  decay  time  constant.  If  another 
presynaptic  spike  arrives  while  a  residual  synaptic  current  remains,  the  conductance  is  not 
allowed  to  increase  above  the  level  assigned  for  a  single  spike  (Ekeberg  et  al.  1991).  The  model 
is  capable  of  producing  action  potentials  with  a  realistic  shape  and  in  a  realistic,  for  the  motor 
cortical  cells,  firing  frequency  range  between  10  and  70  impulses  per  second. 

To  make  the  model  specific  to  a  description  of  the  motor  cortical  network  two  additional 
factors  are  introduced.  First,  each  f  neuron  receives  an  extra  stimulation  by  intracellular  current 
injection.  This  extra  current,  E^,  serves  to  assign  the  preferred  direction  of  the  neuron.  Specific 
expression  for  the  current  is  chosen  in  order  to  reproduce  quantitatively  the  directional  tuning 
curve,  or  the  changes  in  motor  cortical  cell  activity  related  to  changes  of  externally  given 
direction  of  the  movement.  The  experimentally  observed  tuning  curve  can  be  approximated  as  a 
linear  function  of  the  cosine  of  the  angle  between  the  preferred  direction  of  a  cell  and  the 
movement  direction  (Georgopoulos  et  al.  1982;  Schwartz  et  al.  1988).  Within  the  model  this 
condition  leads  to  the  following  relation: 


Ei  =  b  +  acos(^  -a,) 


(1) 


where  the  angle  6  corresponds  to  externally  given  direction,  the  angle  a,-  is  the  preferred  direction 
angle  of  the  i"‘  neuron,  b  and  a  are  adjustable  parameters.  The  observed  tuning  curves  are 
reproduced  if  the  following  values  of  the  parameters  are  used:  6  =  3.0  nA,  a  =  2.2  nA.  Second, 
the  value  of  synaptic  conductance  corresponding  to  the  connection  between /*  (presynaptic)  and 
i"'  (postsynaptic)  neurons  shown  in  Table  1  is  multiplied  by  an  additional  parameter,  the 
coimection  strength  Wy,  and  the  type  of  a  synapse  (excitatory  or  inhibitory)  is  determined  by  the 
sign  of  the  Wy  value  (positive  or  negative  correspondingly).  A  specific  form  of  the  correlation 
between  the  connection  strengths  Wy  and  the  preferred  direction  angles  a,-  is  chosen  in  accordance 
with  the  experimental  findings  (Georgopoulos  et  al.  1993).  It  was  shown  (Georgopoulos  et  al. 
1993)  that  weights  of  functional  connections  between  directionally  tuned  neurons  tend  to  be 
negatively  correlated  with  the  difference  between  the  preferred  directions.  A  first  obvious 
approach  to  model  this  correlation  is  to  represent  an  unknown  function  w,j(a,-a^)  as  the  sum  of 
two  harmonics: 


Wij- 


-■  d  cos(a ,  -  a ;) + c  sin(a ,  -  a  7) 


(2) 


where  d  and  c  are  adjustable  parameters  that  coordinate  relative  contribution  of  symmetric  and 
anti-symmetric  components,  and  d>0.  The  correlation  between  the  weight  of  the  synaptic 
connection  Wy  and  the  angles  a,-,  a,  is  the  main  feature  of  the  suggested  model. 
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Simulations  and  discussion 

The  dynamics  of  the  network  is  governed  by  the  system  of  384  coupled  differential  equations  (8 
equations  per  each  neinon).  Expressions  for  the  extra  currents  and  for  the  synaptic  connection 
strengths  Wy  are  described  by  (2), (3).  The  preferred  direction  angles  «,•  are  uniformly  distributed 
on  the  interval  [-;r,7r].  For  a  given  set  of  parameters  the  system  of  differential  equations  is  solved 
using  fifth-order  Runge-Kutta  formula. 

We  found  that  a  directional  tendency  of  the  network  emerged  due  to  segregation  of 
neurons  imposed  externally  by  extra  currents  F^,  (see  Fig.  1  in  Lukashin  and  Georgopoulos, 

1994).  The  direction  of  the  population  vector  coincided  with  the  externally  given  direction  9  =  0 
deg.  At  the  instant  of  time  150  ms  the  symmetric  component  of  the  synaptic  connections  was 
switched  on,  that  is,  Wy  =  cos(a,-a,).  After  this  moment  the  directional  tendency  becomes  even 
stronger  because  of  an  inhibition  of  the  firing  of  neurons  whose  preferred  direction  are  far  firom 
the  externally  given  direction  0  deg.  The  most  interesting  dynamical  period  begins  at  300  ms, 
when  the  part  of  external  currents  that  realized  the  specification  of  neurons  by  their  preferred 
directions  [second  term  in  the  right-hand  side  of  (1)]  was  switched  off.  Starting  from  this 
moment,  there  is  no  directional  influence  that  would  come  from  outside  the  network,  since  all 
neurons  receive  the  equal  extra  current  =  3.0  nA.  Nevertheless  the  directional  tendency  of  the 
network  is  being  kept  through  the  mechanism  of  the  correlated  interactions  between  cells  (see  the 
spike  trains  and  population  clusters  at  the  interval  300  -  900  ms).  These  simulations  showed  that 
the  experimentally  observed  stability  of  the  directional  tendency  of  the  motor  cortex  lasting 
several  hundred  milliseconds  (Georgopoulos  et  al.  1986;  Georgopoulos  et  al.  1989a; 
Georgopoulos  et  al.  1993)  could  be  ensured  by  specifically  correlated  synaptic  connections 
between  cells  even  without  any  directional  influence  from  areas  outside  the  motor  cortex. 

Finally,  the  ability  of  the  network  to  reproduce  a  basic  type  of  transformations  of  the  neuronal 
population  vector,  namely  its  rotation,  was  also  demonstrated.  The  rotation  of  the  population 
vector  at  90  deg  takes  approximately  100-150  ms  (between  400  and  550  ms),  which  gives  a 
velocity  value  of  about  600-900  deg/s.  This  value  is  quite  close  to  the  experimentally  measured 
angular  velocity  of  the  rotation  of  the  neuronal  population  vector  which  is  about  400-700  deg/s 
(Georgopoulos  et  al.  1989b;  Lurito  et  al.  1991). 

These  results  demonstrate  that  the  main  features  of  experimentally  observed  motor 
cortical  activity  related  to  the  directional  tendency  can  be  reproduced  within  the  suggested 
model.  The  following  remarks  are  noteworthy.  First,  the  results  are  not  specific  to  the  model, 
nor  to  the  set  of  parameters  chosen.  We  checked  the  robustness  of  the  results  in  respect  to 
changes  of  parameters  within  a  reasonable  range,  and  we  obtained  qualitatively  the  same  results 
for  a  model  with  simplified  channel  kinetics  close  to  that  one  suggested  by  (Bush  and  Sejnowski 
1991).  Second,  concerning  the  assignment  of  the  preferred  directions  of  the  neurons,  we  checked 
that  using  external  current  is  not  the  only  way  to  for  this  assignment:  we  obtained  the  same 
results  by  modeling  external  directional  influences  by  adding  extra  neurons  with  correspondingly 
assigned  strengths  of  the  synaptic  connections.  In  the  present  paper  only  one  particular  type  of 
the  correlation  between  the  weight  of  the  synaptic  connection  and  the  difference  between  the 
preferred  directions  of  neuronal  pair  involved  in  the  connection  (3)  has  been  tested.  Previous 
results  of  a  study  of  the  stability  problem  performed  in  (Georgopoulos  et  al.  1993)  within  a  more 
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simple  model  show  that  other  types  of  this  solution  might  be  expected.  A  global  change  of  the 
synaptic  connection  strengths  that  induced  within  the  model  the  rotation  of  the  population  vector 
might  have  its  biological  analog  as  so-called  saliency  systems  which  are  able  to  affect  many 

cortical  synapses  simultaneously  by  release  of  a 
and  Marder  1991;  Tononi  et  al.  1992). 

modulatory  neurotransmitter  (Harris-Warrick 

Table  1.  Parameters  used  for  the  simulations. 

Parameter  Value 

Cell 

Leak  current  equilibrium  potential 

-70  mV 

Sodium  equilibrium  potential 

50  mV 

Potassium  equilibrium  potential 

-90  mV 

Calcium  equilibrium  potential  (voltage  gated) 

150  mV 

Calcium  equilibrium  potential  (NMDA  gated) 

20  mV 

Membrane  capacitance 

0.05  nP 

Membrane  leak  conductance 

0.01  pS 

Sodium  conductance 

1.00  pS 

Potassium  conductance 

0.50  pS 

Calcium-dependent  potassium  conductance 

0.02  pS 

Calcium  accumulation  rate  (voltage  gated) 

0.0040  mV'ms* 

Calcium  accumulation  rate  (NMDA  gated) 

0.0005  mV'ms-* 

Calcium  decay  rate  (voltage  gated) 

0.0300  ms'* 

Calcium  decay  rate  (NMDA  gated) 

0.0030  ms' 

Excitatory  synapses 

Reversal  potential 

OmV 

Conductance  after  spike 

0.003w*  pS 

Decay  time  constant 

100  ms 

Inhibitory  synapses 

Reversal  potential 

-85  mV 

Conductance  after  spike 

0.015W*  pS 

Decay  time  constant 

20  ms 

l.C.  A  model  of  interacting  networks 

Abstract 

The  hypothesis  was  tested  that  learned  movement  trajectories  of  different  shapes  can  be  stored  in, 
and  generated  by,  largely  overlapping  neural  networks.  Indeed,  it  was  possible  to  train  a 
massively  interconnected  neural  network  to  generate  different  shapes  of  internally  stored. 
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dynamically  evolving  movement  trajectories  using  a  general-purpose,  core  part,  common  to  all 
networks,  and  a  special-purpose  part,  specific  for  a  particular  trajectory.  The  weights  of 
connections  between  the  core  units  do  not  carry  any  information  about  trajectories.  The  core 
network  alone  could  generate  externally  instructed  trajectories  but  not  internally  stored  ones,  for 
which  both  the  core  and  the  trajectory-specific  part  were  needed.  All  information  about  the 
movements  is  stored  in  the  weights  of  coimections  between  the  core  part  and  the  specialized 
units  and  between  the  specialized  units  themselves.  Due  to  these  coimections  the  core  part 
reveals  specific  dynamical  behavior  for  a  particular  trajectory  and,  as  the  result,  discriminates 
different  tasks.  The  percentage  of  trajectory-specific  units  needed  to  generate  a  certain  trajectory 
was  small  (2-5%),  and  the  total  output  of  the  network  is  almost  entirely  provided  by  the  core 
part,  whereas  the  role  of  the  small  specialized  parts  is  to  drive  the  dynamical  behavior.  These 
results  suggest  an  efficient  and  effective  mechanism  for  storing  learned  motor  patterns  in,  and 
reproducing  them  by,  overlapping  neural  networks,  and  are  in  accord  with  neurophysiological 
findings  of  trajectory-specific  cells  and  with  neurological  observations  of  loss  of  specific  motor 
skills  in  the  presence  of  otherwise  intact  motor  control . 

Introduction 

Although  a  wealth  of  knowledge  has  accumulated  concerning  the  neural  mechanisms  of  visually 
guided  reaching  (Georgopoulos  et  al.  1993)  and  tracing  (Schwartz  1993)  movements,  and  the 
design  and  performance  of  artificial  neural  networks  for  similar  movements  (see  previous 
sections),  our  knowledge  concerning  the  generation  and  performance  from  memory  of  explicitly 
defined,  learned  movement  trajectories,  such  as  drawing  a  circle,  are  largely  unknown.  Certain 
brain  lesions  can  result  in  apparently  specific  loss  of  particular  motor  skills  ["apraxia",  see  De 
Renzi  1986],  such  as  dressing  or  buttoning  a  garment,  without  affecting  other  motor  skills  (e.g 
driving  a  car)  or  simple  movements  (e.g.  reaching  to  a  target).  It  is  generally  assumed  that 
information  concerning  the  performance  of  the  lost  motor  skill  is  stored  in  the  lesioned  areas 
(commonly  in  the  posterior  parietal  cortex)  or  that  these  areas  are  unique  in  triggering  the 
appropriate  motor  action,  the  motor  pattern  of  which  is  stored  elsewhere.  Whatever  the 
mechanism,  the  crucial  supposition  is  that  the  neural  pattern  of  a  motor  skill  [Lashley’s  "motor 
engram"  ]  is  stored  somewhere  in  toto  so  that,  when  activated,  it  unfolds  in  time  as  a  skilled 
motor  act.  Since  movements  are  the  result  of  interactions  among  neurons  in  the  brain,  it  is 
reasonable  to  hypothesize  that  the  motor  engram  could  be  stored  in  the  set  of  connections  and 
synaptic  strengths  between  interacting  neurons  within  and  among  various  sensorimotor  areas. 
This  distributed  representation  of  the  motor  skill  could  accoimt  for  the  elusiveness  of  the  nature 
and  the  site  of  its  motor  engram.  The  neural  networks  subserving  specific  motor  engrams  could 
be  separate  and  very  specific  in  their  composition  and  coimection  strengths,  so  that  a  particular 
learned,  skilled  action  could  be  generated  by  the  exclusive  activation  of  the  corresponding 
network  with  a  fixed  set  of  cotmection  strengths  (Lukashin  and  Georgopoulos  1994).  A 
generalization  of  this  idea  would  posit  the  existence  of  a  very  large  number  of  specific  networks, 
one  for  each  skill,  a  possible  but  impractical  suggestion.  At  the  other  extreme,  one  and  the  same 
network  could  generate  different  motor  behaviors  by  a  continuous  modulation  of  the  coimections 
in  a  single  network  (Harris- Warrick  and  Marder  1991).  In  this  study  we  entertained  an 
intermediate  hypothesis,  namely  that  learned  motor  skills  are  subserved  by  largely  overlapping 
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networks  with  fixed  connection  strengths.  According  to  this  idea,  the  performance  of  a  learned 
motor  skill  involves  a  network  with  two  kinds  of  units:  (i)  general  purpose,  "core"  units  that  are 
common  to,  and,  therefore,  engaged  with,  all  movements  and  skills,  and  (ii)  very  specialized 
units  that  are  dedicated  to,  and,  therefore,  engaged  with,  only  the  particular  set  of  movement 
trajectories  comprising  a  motor  skill.  Visually  guided  pointing  (1-7)  or  tracing  (8)  movements 
could  be  generated  by  the  core  network,  whereas  learned  skilled  movements  could  be  generated 
by  the  concomitant  activation  of  both  the  core  and  the  specialized  units.  This  would  be  a 
distributed  mechanism  by  which  great  specificity  could  be  achieved  with  a  minimum  of 
dedicated  neural  resources. 

Model 

We  tested  the  hypothesis  above  by  using  massively  interconnected  neural  networks  modeled 
according  to  the  results  of  experimental  studies  (see  Georgopoulos  et  al.  1993);  namely  (i)  the 
units  of  the  network  were  assigned  preferred  directions,  (ii)  the  time-varying,  dynamically 
evolving  outcome  of  the  network  operation  was  calculated  as  the  sum  of  the  vectorial 
contribution  of  these  units  [i.e.  network  population  vector],  and  (iii)  such  population  vectors  were 
added  successively  tip-to-tail  to  create  a  "neural"  trajectory  (Georgopoulos  et  al.  1988). 
Specifically,  if  C,  is  the  unit  preferred  direction  vector  for  the  f  cell,  then  the  neuronal 
population  vector  P  is  defined  as  the  weighted  sum  of  these  vectors: 

P(t)  =  EVi(‘)G  (1) 


where  the  weight  Vff)  is  the  activity  (frequency  of  discharge)  of  the  1“'  unit  at  time  bin  t.  In 
accordance  with  experimental  data  (Georgopoulos  et  al.  1993)  the  preferred  directions  were 
randomly  and  uniformly  distributed  in  space.  A  neural  trajectory  was  obtained  by  attaching 
successive  population  vectors: 

R(fr)  =  EP(t„)  (2) 


where  the  radius-vector  R(ti)  defines  the  point  at  the  neural-vector  trajectory  taken  at  time  bin  4. 

To  calculate  a  neural  trajectory  a  standard  set  of  resistance-capacitance  equations 
governing  the  interactions  between  units  and  their  dynamic  evolution  was  used  (Hopfield  and 
Tank  1986) .  The  time-dependent  output  activity  of  the  f  unit  V,{t)  was  calculated  as  Vft)  = 
tanh[M,(0]j  where  the  variable  m,(0  represents  the  internal  state  (for  instance,  soma  membrane 
potential)  of  the  unit.  The  dynamic  evolution  of  the  pattern  of  activity  of  N  interconnected  units 
was  governed  by  the  following  system  of  equations: 


N 

(t)  +  XwijVj(t)  +  cos(0  -ai) 


j=i 
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where  /  =  1,...,  N;  argument  t  is  shown  for  values  which  depend  on  time;  t  is  a  characteristic  time 
constant;  Wy  is  the  connection  strength  between  units  (/  i).  The  third  term  on  the  right  hand 

side  of  Eq.  3  serves  to  assign  a  preferred  direction  to  the  f  unit  via  extrinsic  input.  The  angle  9 
corresponds  to  an  externally  given  initial  direction,  and  the  angle  a,-  is  regarded  as  the  preferred 
direction  of  the  i"'  unit.  In  the  two-dimensional  case  the  angle  a,-  umquely  defines  the  umt 
preferred  direction  vector  C,-.  In  routine  calculations,  Eq.  3  were  solved  as  an  initial  value 
problem  using  fourth  order  Runge-Kutta  formula  with  automatic  control  of  the  step  size  during 
the  integration. 

The  network  was  trained  to  generate  four  different  two-dimensional  complex  trajectories: 
a  clockwise  circle,  an  orthogonal  bend,  a  counterclockwise  circle  and  a  sinusoid.  Each  of  these 
trajectories  was  generated  by  a  network  that  comprised  a  general-purpose,  core  part  (large  ovoid) 
plus  a  special-purpose  set  of  units  specific  for  the  particular  trajectory  (one  of  the  small  ovoids). 
The  core  part  is  common  to,  and  shared  by,  all  four  networks,  and  is,  therefore,  activated 
regardless  of  the  shape  of  the  trajectory:  the  particular  shape  of  a  trajectory  depends  on  the 
specific  set  of  units  activated,  together  with  the  core  units,  while  the  remaining  trajectory- 
specific  sets  are  inactive.  The  connection  strengths  among  the  units  of  the  core  part  are  fixed  and 
remain  the  same  for  all  trajectories.  On  the  other  hand,  the  connection  strengths  between  the 
core  units  and  the  trajectory-specific  imits,  and  those  among  the  specific  umts  themselves,  are 
allowed  to  change  during  training  {^'variable"  connection  strengths)  of  the  network. 

Training  procedure 

To  train  the  network  to  generate  desired  trajectories  the  variable  synaptic  weights  were  adjusted 
by  means  of  the  simulated  annealing  algorithm  (Kirkpatrick  et  al.  1983).  Specifically,  the 
simulated  annealing  procedure  was  used  to  minimize  the  root-mean-square  (RMS)  error  between 
the  desired  trajectory  shape  and  that  generated  by  the  network: 


F  = 


sl/2 


where  the  radius- vectors  i?/4)  show  the  corresponding  points  at  the  desired  trajectory 

and  at  actual  trajectory  generated  by  the  network  taken  at  time  and  Ra(to)  =  R/tf}.  In  routine 
calculations  we  considered  200  points  {K=  200  in  Eq.4).  Each  step  of  the  simulated  annealing 
procedure  included  a  random  change  of  one  of  the  variable  synaptic  weights  followed  by  an 
entire  recalculation  of  the  trajectory  generated  by  the  network.  The  new  value  of  the  synaptic 
weight  was  accepted  not  only  for  changes  that  lowered  the  RMS  error,  but  also  for  changes  that 
raised  it.  The  probability  of  the  latter  event  was  chosen  such  that  the  system  eventually  obeyed 
the  Boltzmann  distribution  at  a  given  "temperature",  if  the  RMS  error  is  treated  as  the  "energy" 
of  the  system.  The  temperature  was  decreased  according  to  a  cooling  schedule  T„^,  =PT„,  where 
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T„  was  the  temperature  at  the  n"*  step  and  the  value  1  -  yff  was  varied  within  the  interval  from 
5x1 0"^  to  10■^  Each  trial  of  the  training  procedure  was  repeated  with  different  cooling  schedules 
(different  values  of  the  parameter  p)\.o  avoid  the  local  minima  problem.  Generally,  if  the 
cooling  is  sufficiently  slow  for  equilibrium  to  be  established  at  each  temperature,  the  global 
minimum,  i.e.  F  =  0,  can  be  reached  in  the  limit  of  zero  temperature.  We  checked  the  robustness 
of  the  results  with  respect  to  different  series  of  random  numbers  used  during  the  realization  of  the 
simulated  annealing  procedure.  The  amount  of  time  required  to  train  the  network  depended  on 
the  number  of  units  in  the  simulation,  on  the  trajectory  used,  and  on  the  particular  set  of 
connection  weights  among  the  core  units.  For  example,  the  computer  time  required  on  a  single 
YMP  C-90  processor  ranged  from  1.5  to  7.3  CPU  hours  for  each  trial  using  100  units  and  the 
clockwise  circular  trajectory. 

Results  of  simulations  and  discussion 

Given  the  overlapping  design  of  the  networks,  we  sought  to  determine  the  minimal  number  of 
trajectory-specific  units  needed  to  generate  a  particular  trajectory.  Consider  a  core  network 
consisting  ofN^  units.  We  first  assigned  mdjixed  the  synaptic  weights  between  the  units  of  the 
core  part.  These  synaptic  weights  were  assigned  randomly  by  the  relation 


f  2 
1 — arcos 

V 


(CiCj) 


where  jy  was  a  random  number  uniformly  distributed  on  the  interval  [-0.5, 0.5].  The  second  term 
on  the  right  hand  side  of  the  Eq.  5  was  introduced  to  provide  a  negative  correlation  between 
synaptic  weights  and  the  difference  between  preferred  directions  of  the  connected  units  This 
type  of  correlation  between  the  synaptic  weights  and  properties  of  directionally  tuned  units  was 
observed  in  experimental  (Georgopoulos  et  al.  1993)  and  modeling  studies  (Lukashin  and 
Georgopoulos  1994).  Then  a  trajectory-specific  part  with  a  given  number  of  imits  (n  =  1,2,3..) 
was  added  and  the  variable  synaptic  weights  were  adjusted  (see  above)  until  the  network  subset 
(i.e.  core  plus  specific  units)  generated  the  desired  trajectory,  or  imtil  training  failure  was  evident. 
We  considered  the  training  procedure  successful  if  the  resulting  network  was  able  to  generate 
trajectory  which  provided  the  value  of  the  RMS  error  (Eq.  4)  equal  to  or  less  than  unity.  If  this 
was  not  achieved  after  3x10^  steps  of  the  simulated  annealing  procedure,  the  training  was 
considered  unsuccessful.  We  did  not  suppose  any  specificity  of  the  values  of  variable  synaptic 
weights  in  comparison  with  the  weights  of  connections  between  the  core  units.  During  the 
simulated  annealing  procedure  a  new  probe  value  for  variable  synaptic  weight  was  randomly 
selected  in  accordance  with  Eq.  4.  This  means  that  the  probe  values  for  variable  S5maptic  weights 
obeyed  the  same  distribution  function  as  the  fixed  synaptic  weights  for  connections  among  the 
core  units. 
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The  value  of  n  just  sufficient  for  successful  training  was  considered  to  be  the  minimal 
sufficient  number  of  trajectory-specific  units,  n,,  for  this  trial;  this  number  varied  somewhat  from 
trial  to  trial.  For  a  given  number  of  units  in  the  core  part,  the  procedure  above  was  repeated 

ten  times  using  different  but  fixed  synaptic  weights  for  connections  among  the  core  units.  This 
was  carried  out  for  each  type  of  trajectory  and  for  values  ranging  from  25  to  205. 

The  number  of  trajectory-specific  units  sufficient  to  generate  a  particular  trajectory  was 
small:  for  100  core  units,  2-5  trajectory-specific  units  were  sufficient.  It  is  noteworthy  that  this 
finding  is  independent  of  the  particular  shape  of  a  trajectory.  Moreover,  the  ratio  of  trajectory- 
specific  units  over  the  number  of  core  units  {n/N^  decreased  as  the  number  of  core  units 
increased.  Although  these  results  cannot  be  directly  extrapolated  to  very  large  networks,  larger 
simulations  could  yield  either  no  improvement  in  accuracy  (saturation)  or  further  decreases  in 
trajectory  error.  In  the  latter  case  the  fraction  of  the  trajectory-specific  units,  relative  to  the 
number  of  core  units,  could  be  even  less  than  our  estimate  of  2-5%. 

Thus,  in  the  framework  of  our  model  the  larger  core  part  of  the  network  does  not  carry 
any  information  about  possible  movements  in  the  static  state  because  the  weights  of  connections 
between  the  core  units  are  the  same  for  all  trajectories.  All  information  about  the  movements  is 
stored  in  the  weights  of  connections  between  the  core  part  and  the  specialized  miits,  and  between 
specialized  units  themselves.  However,  once  the  dynamics  gets  started  by  activation  of  one  of  the 
specialized  set,  the  core  part  reveals  specific  dynamical  behavior  for  a  particular  trajectory,  due 
to  the  driven  forces  from  the  specialized  units.  Therefore,  during  the  dynamics  the  core  part  does 
discriminate  different  tasks.  Note,  that  the  core  part  also  actively  influences  the  dynamics 
through  the  feedback  connections  to  the  specialized  units.  The  roles  played  by  the  core  part  and 
by  specialized  parts  are  the  following.  Since  the  number  of  specialized  units  are  negligibly  small 
in  comparison  with  the  size  of  the  core  part,  the  total  output  of  the  whole  network  is  almost 
entirely  provided  by  the  dynamical  behavior  of  the  core  network  which  can  translate  the 
information  to  the  lower  levels  of  the  central  nervous  system.  The  role  of  the  small  specialized 
parts  is  to  receive  information  about  the  beginning  of  the  movement  and  to  drive  the  dynamical 
behavior  of  the  whole  network. 

Recent  neurophysiological  studies  (e.g.  Ashe  et  al.  1994)  have  shown  that  a  small  percent 
(1-10%)  of  cells  recorded  during  performance  of  learned  movements  from  memory  are  very 
specific  to  a  particular  trajectory,  whereas  a  relatively  large  number  of  cells  are  engaged  during 
both  simple  pointing  movements  and  diuing  performance  of  the  specialized  movements.  These 
observations  are  in  close  quantitative  agreement  with  the  results  of  the  present  study;  indeed,  the 
experimental  results  can  be  regarded  as  reflecting  the  limit  of  the  theoretical  results  obtained  in 
this  study.  In  the  brain,  specialized  networks  could  be  activated  by  various  cortical  and 
subcortical  structures  including  the  cerebellum  and  basal  ganglia.  Finally,  the  architectiue  of 
overlapping,  massively  interconnected  networks  with  a  minimum  of  specialized  units  could  be 
useful  to  other  applications  requiring  the  production  of  very  specific  outcomes:  this  architecture 
is  efficient  and  effective,  for  it  maximizes  the  specificity  that  can  be  obtained  while  minimizing 
the  number  of  specific  units  and  allowing  for  a  common  core  to  be  shared  by  different 
applications. 
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2.  Development  of  a  model  of  interactions  between  the  cortical  and  the  spinal  networks 
dealing  with  generating  time-varying  motoneuronal  outputs  for  movements  in  space 
(publications  7-8). 

These  studies  comprised  (A)  development  of  a  “spinal-type”  network  processing  force 
production,  and  (B)  modeling  of  the  interactions  between  motor  cortical  and  spinal  networks. 

2.A.  A  spinal-type  network  processing  force  production 


Abstract 

We  have  developed  a  model  that  simulates  possible  mechanisms  by  which  supraspinal  neuronal 
signals  coding  forces  could  converge  in  the  spinal  cord  and  provide  an  ongoing  integrated  signal 
to  the  motoneuronal  pools  whose  activation  results  in  the  exertion  of  force.  The  model  consists 
of  a  three-layered  neural  network  coimected  to  a  two-joint-six-muscle  model  of  the  arm.  The 
network  layers  represent  supraspinal  populations,  spinal  cord  intemeurons,  and  motoneiuonal 
pools.  We  propose  an  approach  to  train  the  network  so  that,  after  the  synaptic  coimections 
between  the  layers  are  adjusted,  the  performance  of  the  model  is  consistent  with  experimental 
data  obtained  on  different  organisms  using  different  experimental  paradigms:  the  stifftiess 
characteristics  of  human  arm;  the  structure  of  force  fields  generated  by  the  stimulation  of  the 
fi'og's  spinal  cord;  and  a  correlation  between  motor  cortical  activity  and  force  exerted  by  monkey 
against  an  immovable  object.  The  model  predicts  a  specific  pattern  of  connections  between 
supraspinal  populations  coding  forces  and  spinal  cord  intemeurons:  the  weight  of  coimection 
should  be  correlated  with  directional  preference  of  interconnected  units.  Finally,  our  simulations 
demonstrate  that  the  force  generated  by  the  sum  of  neural  signals  can  be  nearly  equal  to  the 
vector  sum  of  forces  generated  by  each  signal  independently,  in  spite  of  the  complex 
nonlinearities  intervening  between  supraspinal  commands  and  forces  exerted  by  the  arm  in 
response  to  these  commands. 


Introduction 

Recordings  of  neuronal  activity  in  the  motor  cortex  of  behaving  animals  have  shown  that  the 
direction  and  overall  trajectory  of  arm  movements  are  markedly  related  to  the  activity  of  motor 
cortical  cells  (Georgopoulos  et  al.  1982;  Georgopoulos  et  al.  1986;  Kalaska  et  al.  1989;  Caminiti 
et  al.  1990;  Hocherman  and  Wise  1991;  Georgopoulos  et  al.  1993;  Schwartz  1994).  These 
observations  have  been  extended  to  the  study  of  forces  exerted  by  the  arm  against  an  immovable 
object  (Georgopoulos  et  al.  1992).  A  monkey  was  trained  to  exert  forces  on  an  isometric  handle 
in  the  presence  of  a  constant  force  bias.  First  the  monkey  was  required  to  exert  a  postural  (static) 
force  P,  which  compensated  a  given  bias  force  B.  After  a  holding  period,  a  cue  instructed  the 
monkey  to  exert  a  force  S  so  that  the  net  force  acting  on  the  handle  (i.e.,  the  force  exerted  by  the 
monkey  S  plus  the  bias  force  applied  to  the  handle  B)  was  in  a  visually  specified  direction  N. 
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Note,  that  the  net  force  is  equivalent  to  the  incremental  (dynamic)  component  I  of  the  force 
exerted  by  the  subject:  I=SBP=S+B=N .  Eight  net  force  directions  and  eight  bias  force  directions 
were  employed.  Two  principal  findings  were  reported;  (i)  the  activity  of  single  cells  showed 
approximately  the  same  directional  tuning  properties  when  the  arm  exerts  a  force  without 
moving  as  when  it  moves  through  space  (Georgopoulos  et  al.  1993),  and  (ii)  as  a  population,  the 
activity  of  cells  reflected  the  direction  of  the  incremental  force  I  and  not  the  actual  force  S 
exerted  by  the  subject.  Accordingly,  it  was  hypothesized  (Georgopoulos  et  al.  1992; 
Georgopoulos  1994)  that  the  postural  and  incremental  signals  are  controlled  by  different  neural 
systems,  and  that  these  signals  would  converge  in  the  spinal  cord  and  produce  a  resulting 
integrated  signal  to  the  motoneuronal  pools. 

The  hypothesis  concerning  an  integration  of  separate  force  signals  imphes  a  "linear 
summation  rule":  the  postural  signal,  which  generates  the  force  P,  and  the  incremental  signal, 
which  generates  the  force  I,  should  converge  at  the  level  of  spinal  cord  in  such  a  way  that  the 
integrated  signal  would  generate  the  force  S  that  is  equal  to  the  vector  sum  P+I.  Only  in  this  case 
the  net  force  acting  on  the  handle  {S+B)  will  be  equal  to  the  desired  net  force  N.  Although  an 
appropriate  summation  of  neuronal  signals  is  conceivable  (Redish  and  Touretzky  1994),  the 
linear  summation  of forces  seems  unlikely,  due  to  the  complex  nonlinearities  that  characterize 
the  mechanical  properties  of  limbs,  and  the  interactions  both  among  neurons  and  between 
neurons  and  muscles  (Bizzi  et  al.  1991;  Kalaska  and  Crammond  1992).  However,  recent 
evidence  derived  by  focal  microstimulation  of  the  fi*og's  spinal  cord  (Bizzi  et  al.  1991;  Giszter  et 
al.  1993;  Mussa-Ivaldi  et  al.  1994)  has  revealed  that  the  simultaneous  activation  of  two  distinct 
spinal  sites  leads  to  the  vectorial  summation  of  the  end-point  forces  generated  by  the  stimulation 
of  each  site  separately.  Although  spinal  microstimulation  studies  have  not  been  performed  in 
primates,  it  is  reasonable  to  hypothesize  a  similar  plan  of  spinal  organization.  Then  the  synaptic 
coimections  between  the  supraspinal  populations  and  the  spinal  cord  could  provide  a  concomitant 
activation  of  a  number  of  the  spinal  sites  associated  with  different  force  fields,  which  are 
summed  in  a  linear  fashion. 

The  purpose  of  the  present  article  is  to  propose  and  analyze  a  model  that  suggests 
possible  mechanisms  by  which  the  supraspinal  commands  could  be  integrated  at  the  spinal  cord 
level  and  translated  into  exertion  of  a  required  force.  The  model  consists  of  a  three-layered  neural 
network  connected  to  a  two-joint-six-muscle  model  of  the  arm.  The  layers  of  the  network 
represent  supraspinal  neuronal  populations,  spinal  cord  mtemeurons,  and  motoneuronal  pools, 
respectively.  The  key  idea  is  to  train  the  network  so  that  the  model  reproduces  quantitatively 
experimental  data  for  stiffness  characteristics  of  human  arm  (Mussa-Ivaldi  et  al.  1985;  Shadmehr 
et  al.  1993;  Tsuji  et  al.  1995),  and  then  to  test  the  model  by  comparing  its  performance  with 
experimental  data  obtained  from  recordings  of  cell  activity  in  the  motor  cortex  of  monkey 
(Georgopoulos  et  al.  1992)  and  by  microstimulations  of  the  fi'og's  spinal  cord  (Bizzi  et  al.  1991; 
Giszter  et  al.  1993;  Mussa-Ivaldi  et  al.  1994).  In  particular,  the  model  reconciles  the 
nonlinearities  that  intervene  between  supraspinal  commands  and  forces  exerted  by  the  arm  with 
the  hypothesis  concerning  the  linear  summation  of  forces  encoded  in  activities  of  potentially 
separated  supraspinal  neuronal  ensembles. 
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Model 


Neural  network  model.  In  the  model,  the  direction  and  magnitude  of  incremental  force  / 
is  encoded  in  activities  of  directionally  tuned  motor  cortical  cells  (Georgopoulos  et  al.  1992). 
Although  the  supraspinal  representation  of  postural  force  P  is  less  known  (see,  however,  Kalaska 
et  al.  1989;  Wise  1993;  Georgopoulos,  1994),  we  assume  the  same  distributed  representation  as 
for  the  incremental  force.  Each  population  consists  of  n  supraspinal  (SS)  units  of  m  neurons  each. 
All  neurons  belonging  to  the  same  SS  imit  possess  the  same  preferred  directions  C,-  (i=l,...,n) 
with  cosine-like  timing  function  (see  Equation  1  below),  and  the  preferred  directions  C,-  are 
uniformly  distributed  in  2-D  space  (Georgopoulos  et  al.  1993;  Georgopoulos  1994).  Therefore,  if 
the  activity  of  supraspinal  population  represents  a  force  signal  F  (F=I  or  P),  then  the  activity  of 
the  z'-th  SS  unit  is  given  by  the  expression: 

)=Y(i+cosepcJ 


where  coefficient  Up  is  proportional  to  the  number  of  cells  recruited  to  represent  the  magnitude  of 
force  F  (i.e.,  the  larger  is  the  magnitude  of  the  encoded  force,  the  more  neurons  from  a  given  SS 
unit  are  recruited  (Evarts  et  al.  1983;  Redish  and  Touretzky,  1994)),  and  6pc  is  the  angle  between 
the  direction  of  force  F  and  the  unit's  preferred  direction  C. 

All  units  from  both  supraspinal  populations  are  connected  to  the  four  units  at  the  IN  layer 
representing  the  level  of  spinal  cord  intemeurons.  To  calculate  activities  of  the  iW  units  VN 
(/=1,...,4)  we  use  the  following  sigmoid  activation  function: 

F/''  =  |t  +  tanh(Tj  +  Uf  )] 


where  7)  is  a  synaptic  input  furnished  by  sources  other  than  the  populations  coding  forces.  These 
signals  provide  an  initial  pattern  of  activity,  needed  to  establish  an  initial  position  of  the  arm  (see 
Section  3  below),  before  the  commands  coding  forces  P  and  I  reach  the  IN  layer.  The  synaptic 
input  Uj^  generated  by  supraspinal  force  signals  is  given  by  the  sum: 

Uf  =  E  Wji  [Pf  (  P )  +  Vf  ( I )] 


where  w,,  is  the  weight  of  connection  between  the  z-th  unit  from  the  population  coding  force  P  or 
/  and  the y-th  unit  belonging  to  the  IN  layer  (for  the  sake  of  simplicity  we  assume  the  same 
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connectivity  matrix  for  both  supraspinal  populations). 

Each  of  the  IN  units  is  connected  to  all  six  units  at  the  MV  layer  representing  the  level  of 
motoneuronal  pools.  Activities  of  A/A’ units  (A=l,...,6)  are  calculated  using  the  relation: 

F/"'  =  ^[i  +  tanh(ur)] 


where  is  a  synaptic  input  generated  by  /A  units: 

4 

f.IN  —  X'  jrIN 

Uk  -Z^ZkjVj 


where  is  the  weight  of  connection  between  the y-th  unit  from  the  IN  layer  and  the  A:-th  unit 
from  the  MN  layer. 

Model  of  the  arm.  We  use  a  planar  two-joint-six-muscle  model  utilized  in  (Hogan  1985; 
Flash  1987;  Katayama  and  Kawato  1993).  The  arm  is  modeled  by  two  rigid  segments  ("upper 
arm  "  and  "forearm")  of  equal  length.  The  upper  arm  is  attached  proximally  to  an  immovable 
"clavicle"  via  the  shoulder  joint  and  distally  to  the  forearm  via  the  elbow  joint.  Both  single-  and 
double-joint  muscles  are  included.  One  end  of  both  flexor  and  extensor  single-joint  muscles 
controlling  the  shoulder  joint  is  attached  to  the  clavicle  at  a  distance  b  from  the  shoulder  joint 
(flexor  and  extensor  are  attached  at  different  sides  of  the  joint).  The  other  ends  of  these  muscles 
are  attached  to  the  upper  arm  segment  at  a  distance  L  from  the  shoulder  joint.  Single-joint 
muscles  controlling  the  elbow  joint  are  attached  to  the  upper  arm  segment  at  the  distances  b  from 
the  elbow  joint  and  to  the  forearm  segment  at  a  distanced  from  the  elbow  joint.  Two-joint 
muscles  are  attached  at  distances  b  from  each  joint.  Note,  that  the  muscle  attachment  implies  that 
the  moment  arms  of  muscles  depend  upon  the  joint  angles.  Below  we  assume  that  the  distance 
between  the  shoulder  and  elbow  joints  is  equal  to  L,  and  use  the  limit  b«L\  specifically,  we  set 
6=0.01  m,  L=0.33  m. 

All  muscles  are  modeled  by  nonlinear  springs  with  the  following  length-tension 
relationship  (Feldman  1966;  Shadmehr  and  Arbib  1992): 

/^=aj^exp[p  dk-lOl-l}  (lk>in 

fk  =  0  (ik^lO 


where is  a  contraction  force  developed  by  the  ^-th  muscle,  f  is  the  actual  muscle  length  and  If 
is  the  muscle  intrinsic  rest  length;  parameters  a  and  f  are  constants,  and  a=10  N, P=\00  m'*. 
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These  specific  values  of  parameters  a  and  P  were  chosen  to  fit  by  equation  (6)  the  experimental 
data  of  Feldman  (1966,  as  presented  in  Fig.3C  of  the  paper  by  Shadmehr  and  Arbib  1992). 

Interactions  between  neurons  and  muscles.  Each  unit  firom  the  MN  layer  innervates  one 
muscle  and  controls  the  muscle  intrinsic  rest  length  4°  (Feldman  1966;  Shadmehr  and  Arbib 
1992).  We  use  the  following  relation  between  4“  and  activity  of  motoneuronal  unit: 

^  0  _  I  max  _j_  'y  MN  ^  j  min  _  ^  max  ^ 


where  and  are  the  minimal  and  maximal  muscle  rest  length,  respectively,  and 
pm_0  26  m,  /'"'“=0.30  m.  Equation  (7)  implies  that  the  higher  is  the  activity  V^,  the  shorter  is 
the  muscle  rest  length  4°- 

Equilibrium  position  and  restoring  force  field.  Given  a  particular  set  of  muscle  rest 
lengths  If  {h=\,. a  corresponding  equilibrium  configuration  of  the  arm  can  be  calculated  as 
follows.  Let  a  configuration  of  the  arm  be  defined  by  a  pair  of  shoulder  and  elbow  joint  angles 
((Pj  and  (pg,  respectively;  0°<(pj<135“,  0°<(p^l80°).  Then  the  set  of  actual  muscle  lengths  4 
(^1,...6)  is  uniquely  defined  by  the  geometry  of  muscles  attachment.  If  muscle  rest  lengths  If 
and  muscle  actual  lengths  4  are  known,  the  forces  produced  by  muscles  can  be  calculated  using 
(6).  Let  (/y  Y  be  the  vector  of  muscle  forces  (subscripts  1  and  2  refer  to  the  single-joint 

shoulder  flexor  and  extensor,  subscripts  3  and  4  refer  to  the  elbow  single  joint  flexor  and 
extensor,  and  subscripts  5  and  6  refer  to  the  double-joint  flexor  and  extensor).  The  transformation 
between  the  vector  of  muscle  forces  and  the  vector  of  net  joint  torques  (t^  4  )^»  where  4  and  4  are 
shoulder  and  elbow  torques,  respectively,  is  given  by  the  matrix  of  muscle  moment  arms.  For  the 
model  described  above,  this  relationship  is: 
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where  r=bsm^)^,  r=bsm(^^.  The  equilibrium  configuration  of  the  arm  is  defined  by  the  condition 
that  both  net  joint  torques,  and  4,  are  equal  to  zero.  Consequently,  the  equilibriiim  joint  angles 
(p,®*  and  (p/^  can  be  obtained  as  the  solution  of  the  system  of  two  nonlinear  equations  (8)  under 
conditions  4=0  and  4=0. 

The  force  exerted  by  the  end-point  of  the  arm  is  the  result  of  difference  between  current 
and  equilibrium  end-point  positions.  Generally,  each  equihbrium  position  is  characterized  by  a 
field  of  restoring  forces  generated  by  muscles  in  response  to  displacements  of  end-point  firom 
equilibrium  position.  For  any  given  set  of  muscle  rest  lengths,  the  restoring  force  field  can  be 
calculated  as  follows.  Let  a  current  end-point  position  be  defined  by  joint  angles  and  7^.  Then 

Equation  (8)  can  be  used  to  calculate  the  joint  torques  4  and  4  (in  particular,  if  the  arm 
configuration  is  such  that  9^=9/’  and  9^=9/’,  Equation  8  gives  t=0  and  4=0).  The 
transformation  between  net  joint  torques  and  x-y  components  of  restoring  end-point  force  (Qj  and 
Qy)  is  given  by  the  inverse  Jacobian  matrix: 


I 

K  ^yj 
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The  vector  field  {S^(9j,9e),  Qy(9i,9e)}  is  the  restoring  force  field  associated  with  the  given  set  of 
muscle  rest  lengths. 

Performance  of  the  model.  All  parameters  related  to  the  arm  construction  and  elastic  properties 
of  muscles  will  be  fixed  throughout  the  rest  part  of  the  paper.  Since  the  restoring  force  field  and 
equilibrium  end-point  position  are  uniquely  determined  by  activities  of  the  MV  units  (6-9),  the 
relation  between  force  signals  encoded  in  supraspinal  populations  (input)  and  the  force  exerted 
by  the  arm  (output)  is  determined  by  synaptic  connections  and  z^j  between  neurons  in  different 
layers  (Equations  1-5).  Therefore,  a  desired  performance  of  the  model  can  be  attained  by  training 
of  the  network,  i.e.,  by  an  adjustment  of  synaptic  connections. 

Training-testing  procedure 

In  the  present  section  of  the  paper  we  describe  in  detail  the  method  used  to  adjust  the  synaptic 
connections.  First  the  synaptic  connections  z^j  between  the  IN  and  MN  layers  have  been  adjusted 
so  that  any  pattern  of  activity  at  the  IN  layer  produces  a  restoring  force  field  whose  structure  is 
consistent  with  experimental  data  obtained  for  human  arm.  In  particular,  each  /Aimit  is 
attributed  to  a  specific  point  in  the  workspace,  which  is  the  equilibrium  end-point  position 
produced  by  the  activation  of  this  unit  ("directional  preference"  of  the  unit  given  a  reference 
position  in  the  workspace).  Then  the  synaptic  connections  Wy,  between  supraspinal  populations 
and  the  Z/V  units  have  been  assigned  so  that  activities  of  the  SS  units  weighted  in  accordance  with 
their  preferred  directions  C,-  are  transformed  into  activities  of  the  /?V  units  weighted  in  accordance 
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with  their  directional  preferences  in  the  workspace  of  the  arm  (connections  between  the  /AT  and 
MNvants  are  kept  fixed  to  retain  the  structure  of  force  fields).  Finally,  the  overall  performance  of 
the  model  have  been  tested  by  independent  or  simultaneous  activations  of  supraspinal 
populations. 

Connections  between  the  IN  and  MN  units.  Forces  exerted  by  individual  muscles  depend 
upon  the  muscle  rest  lengths  (6),  which  are  controlled  by  activities  of  the  corresponding  MN 
units  (7).  Therefore,  activities  of  different  MA^  units  should  be  correlated  to  produce  muscle 
synergies  that  ensure  desired  characteristics  of  restoring  force  fields.  Since  any  unit  firom  the  IN 
layer  is  connected  to  all  units  at  the  MN  layer  (4,5),  correlated  patterns  of  the  MN  activity  can  be 
provided  by  a  set  of  "synergetic"  synaptic  connections  between  the  IN  and  A/N  units.  The 
problem  we  address  is  to  find  a  set  of  connections  z,^j  such  that  the  restoring  force  fields  generated 
by  the  model  would  be  similar  to  the  fields  measured  for  human  arm  (Mussa-Ivaldi  et  al.  1985; 
Shadmehr  et  al.  1993;  Tsuji  et  al.  1995).  In  these  studies,  human  arm  stiffiiess  characteristics 
have  been  determined  by  measuring  the  restoring  forces  for  displacements  firom  an  equilibrium 
end-point  position.  The  stiffiiess  was  represented  graphically  by  an  ellipse,  characterized  by  its 
size,  shape,  and  orientation.  The  results  indicated  that  the  shape  and  orientation  of  the  stiffiiess 
ellipse  are  strongly  dependent  on  arm  configuration  in  an  orderly  fashion  (these  parameters, 
however,  remain  invariant  among  subjects  and  over  time).  To  adjust  the  performance  of  the 
model  to  these  experimental  data  we  used  the  following  procedure. 

Four  equilibrium  positions  of  the  arm  were  chosen  in  the  workspace.  For  each  position, 
we  calculated  the  set  of  muscle  rest  lengths  4°  (A=l,...,6)  ensuring  the  size,  shape  and  orientation 
of  the  stiffiiess  ellipse  in  close  agreement  with  these  characteristics  measured  experimentally 
(Mussa-Ivaldi  et  al.  1985)  for  the  human  arm  at  a  similar  location  in  the  workspace.  Note,  that 
due  to  the  presence  of  more  muscles  than  there  are  joints,  different  sets  of  muscle  rest  lengths 
may  result  in  the  same  equilibrium  position  of  the  arm.  However,  it  is  this  redundancy  that 
makes  it  possible  to  find  a  required  set  of  muscle  rest  lengths  imposing  the  stiffiiess 
characteristics  as  additional  conditions.  Starting  fi"om  (8)  we  have  derived  explicit  relations 
between  muscle  rest  lengths  and  the  stiffhess  characteristics  of  the  arm  using  straightforward 
algebra,  cumbersome  as  it  is  (Appendix  A). 

Each  of  these  equilibrium  positions  was  attributed  to  one  of  the  EV  units:  namely,  for  each 
of  four  IN  units  (/=1,...,4),  six  parameters  z^j  (A=l,...,6)  were  adjusted  so  that  the  maximal  activity 
of  the y-th  /A  unit  (whereas  other  units  are  inactive)  generates  the  attributed  equilibrium  position. 

The  resulting  connectivity  matrix  z^y  (A=l,...,6;y=l,...,4)  was  fixed  and  used  during  the 
following  testing  procedure.  Different  patterns  of  IN  activity  were  sampled  to  produce  new 
equilibrium  positions,  for  which  stiffiiess  ellipses  were  calculated.  We  found  that  the  shape  and 
orientation  of  ellipses  strongly  depend  on  arm  configiuation,  reproducing  well  two  main  features 
of  the  human  arm  stiffhess  (Mussa-Ivaldi  et  al.  1985;  Tsuji  et  al.,1995):  (i)  the  direction  of 
maximum  stiffness  at  any  location  is  approximately  oriented  along  a  radial  line  joining  the  hand 
to  the  shoulder,  and  (ii)  the  stiffiiess  ellipse  becomes  more  anisotropic  as  the  hand  approaches 
distal  positions. 

Since  the  stiffiiess  describes  the  relation  between  force  and  displacement  vectors  only  in 
the  vicinity  of  equilibrium,  the  performance  of  the  model  was  further  tested  for  relatively  large 
displacements  from  equilibrium  positions.  The  structure  of  restoring  force  fields  is  similar  to  that 
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measured  experimentally  (Shadmehr  et  al.  1993)  for  the  same  range  of  displacements,  including 
a  nonlinearity  with  respect  to  displacement. 

In  summary,  the  use  of  experimental  data  for  few  equilibrium  positions  as  training 
"examples"  makes  it  possible  to  find  the  connectivity  matrix  z  that  provides  the  required  synergy 
of  muscles  at  other  locations  in  the  workspace. 

Connections  between  the  SS  and  IN  units.  The  connectivity  matrix  w  transforms  a  force 
signal  encoded  in  activities  of  SS  units,  which  are  characterized  by  preferred  directions  (1), 
into  a  signal  encoded  in  activities  of  units,  which  are  attributed  to  specific  equilibrium 
positions  in  the  workspace.  Consider  a  point  chosen  approximately  at  the  center  of  the 
workspace.  Let  Dj  be  a  unit  vector  pointed  fi'om  the  center  of  the  workspace  to  the  equilibrium 
end-point  position  attributed  to  the y-th  ///unit  (in  a  sense,  Dj  is  the  preferred  direction  vector  of 
the y'-th  777  unit).  We  propose  the  following  expression  for  the  connectivity  matrix: 

4 

Wji 


where  is  the  angle  between  vectors  D  and  C,  and  4/n  is  a  normalizing  coefficient  (n  is  the 
number  of  SS  units).  The  idea  underl5dng  the  use  of  the  above  relation  is  that  Equation  10 
provides  a  negative  correlation  between  connection  strength  Wj^  and  angle  between  vectors  Dj  and 
Cj,  i.e.,  the  connection  strength  is  correlated  with  directional  preference  of  connected  units.  The 
cosine  function  (10)  and  the  normalizing  coefficient  are  chosen  to  obtain  activities  ofTA^imits 
Vj^  (2)  in  a  simple  explicit  form.  Indeed,  substituting  (1,3,10)  into  (2),  one  has 

Vj'^  =  ^[j  +  t^r)h{Tj  +  ap  cosOp^.  +  ai  cos0id.  )] 


(to  derive  (1 1)  we  have  used  the  uniformness  of  distribution  of  preferred  directions  C,-  (/=!,...,«) 
and  straightforward  trigonometric  relationships). 

In  summary,  the  cormectivity  matrix  w  (10)  transforms  supraspinal  signals  coding  forces 
P  and  7(1)  into  a  pattern  of  activity  of  the  7A^  units  (1 1)  in  a  way  that  the  closer  is  the  direction  of 
force  to  be  exerted  (either  P  or  7)  to  the  direction  of  vector  Dj,  the  higher  is  the  activity  of  the y-th 
72V  unit. 

Transformation  of  neuronal  signal  into  exertion  of  force 

After  the  synaptic  cormections  Wjf  and  were  adjusted  (see  Section  3)  and  fixed,  we  analyzed 
the  performance  of  the  model  as  follows. 

First  coefficients  ap  and  aj  characterizing  the  magnitude  of  force  signals  were  set  to  zero. 
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and  additional  inputs  T)  (1 1)  were  adjusted  to  obtain  such  initial  pattern  of  the  /iV  activity  that 
produces  the  equilibrium  end-point  position  located  approximately  at  the  center  of  the 
workspace.  Using  this  point  as  the  origin,  the  preferred  direction  vectors  Dj  (/=1,...,4)  were 
defined  as  unit  vectors  pointed  firom  the  origin  to  one  of  four  equilibrium  positions. 

Eight  signals  coding  the  postural  forces  P  of  equal  magnitude  (a,,=0.3)  and  different 
directions  ranging  firom  0  to  360°  were  given,  and  corresponding  patterns  of  the  IN  activity  were 
calculated  in  the  absence  of  incremental  force  signal  (a/=0).  Each  of  these  patterns  produces  a 
new  equilibrium  end-point  position,  which  differs  fi’om  the  initial  position  of  the  arm.  Due  to  this 
difference,  a  restoring  force  appears.  If  the  arm  is  free  to  move,  then  the  restoring  force  would 
result  in  a  new  position  of  the  arm.  Otherwise,  the  restoring  force  can  be  referred  to  the  force 
exerted  by  the  arm  against  an  immovable  object  measured  in  the  experiment  (Georgopoulos  et  al. 
1992).  We  have  calculated  end-point  forces  generated  by  the  eight  neuronal  signals  using  (1-11). 
We  found  that  the  model  properly  transforms  neuronal  signals  of  different  directions  into 
exertion  of  force.  A  difference  between  directions  of  force  signals  and  generated  forces  can  be 
reduced  by  an  additional  adjustment  of  the  connectivity  matrix  w. 

To  test  whether  the  "vectorial  summation  rule"  (see  above)  is  valid  in  the  fi'amework  of 
the  model  we  calculated  restoring  end-point  forces  produced  in  response  to  simultaneous 
activation  of  both  supraspinal  populations  (i.e.,  ap=0.3,  a/=0.3).  We  found  that  in  all  64  cases 
the  sum  of  neuronal  signals  is  transformed  into  the  force  that  is  nearly  equal  to  the  vector  sum  of 
forces  generated  by  each  signal  separately. 

The  performance  of  the  model  can  be  interpreted  in  the  light  of  experimental  data 
reported  in  (Georgopoulos  et  al.  1992;  see  also  Section  1  of  the  present  paper  for  a  discussion). 

To  this  end,  we  refer  the  net  force  N  acting  on  an  immovable  handle  in  the  presence  of  bias  force 
B  (Georgopoulos  et  al.  1992)  to  the  force  calculated  within  the  model  as  the  difference  iSBP, 
where  S  is  the  force  generated  by  the  sum  of  postural  and  incremental  signals,  and  P  is  the  force 
generated  in  response  to  the  postural  signal  (we  suppose  that  the  postural  signal  P  compensates 
exactly  the  bias  force  B,  i.e.,  P=BP).  The  performance  of  the  model  would  be  consistent  with  the 
experimental  findings  (Georgopoulos  et  al.  1992)  if  the  direction  of  "net  force"  5BP  is  close  to 
the  direction  of  incremental  component  /,  regardless  of  the  direction  of  postural  force  P.  In  other 
words,  the  vectorial  summation  rule  S=P+I  should  be  valid.  We  have  calculated  the  net  forces 
SBP  for  eight  different  directions  of  incremental  signal  and  eight  different  directions  of  postural 
signal.  We  found  that  while  the  direction  of  postural  force  P  changes  firom  0  to  360°  within  each 
cluster,  the  direction  of  net  force  iSBP  remains  almost  invariant  and  close  to  the  direction  of 
incremental  component  I.  In  fact,  this  is  another  interpretation  of  the  vectorial  summation  rule. 

Summation  of  restoring  force  fields 

Although  parameters  of  the  model  have  been  adjusted  to  reproduce  the  results  of  specific 
experiments  performed  on  human  and  monkey  (see  above),  it  is  useful  to  compare  the  restoring 
force  fields  generated  by  the  model  with  experiments  performed  on  the  spinalized  firog  (Bizzi  et 
al.  1991;  Giszter  et  al.  1993;  Mussa-Ivaldi  et  al.  1994).  The  results  obtained  by  microstimulation 
of  the  frog's  spinal  cord  suggested  that  the  neural  circuits  in  the  spinal  cord  are  organized  in  a  set 
of  control  modules  that  "store"  a  few  limb  postures  in  the  form  of  convergent  force  fields  acting 
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on  the  limb's  end-point.  Moreover,  it  was  demonstrated  (Bizzi  et  al.  1991;  Mussa-Ivaldi  et  al. 
1994)  that  simultaneous  stimulation  of  two  distinct  spinal  sites  results  in  a  field  of  forces 
proportional  to  the  vector  sum  of  the  fields  induced  by  the  independent  stimulation  of  each  site. 
These  findings  have  led  to  a  new  concept  of  the  motor  control  (Bizzi  et  al.  1991;  Mussa-Ivaldi 
and  Giszter  1992)  based  on  the  vector  combination  of  a  few  convergent  fields  to  produce  a  vast 
repertoire  of  motor  behavior.  Below  we  demonstrate  that  the  proposed  model  generates  restoring 
force  fields  whose  properties  are  in  close  similarity  with  those  reported  in  (Bizzi  et  al.  1991; 
Giszter  et  al.  1993;  Mussa-Ivaldi  et  al.  1994). 

Justification  of  the  model.  To  establish  a  correspondence  between  the  model  and 
experiment  (Bizzi  et  al.  1991;  Giszter  et  al.  1993;  Mussa-Ivaldi  et  al.  1994)  we  used  the 
following  approach.  The  IN  layer  of  the  model  was  discoimected  fi'om  the  supraspinal  layers,  i.e. 
the  model  was  "spinalized".  First  of  all,  we  calculated  the  force  field  generated  under  condition 
that  all  units  at  the  IN  layer  are  inactive.  Since  the  M?/ units  are  also  inactive,  the  muscle  rest 
lengths  in  this  case  are  equal  to  their  maximal  possible  values  (7).  We  refer  this  field  to  the 
"resting"  force  field  measured  in  the  absence  of  stimulation  (Giszter  et  al.  1993).  The  stimulation 
of  a  particular  spinal  site  in  the  experiment  was  referred  to  the  activation  of  a  particular  imit  at 
the  IN  layer.  The  activation  of  distinct  units  resulted  in  different  force  fields  converging  to 
different  equilibrium  points.  We  referred  these  fields  to  the  "total"  force  fields  measured  during 
stimulations  of  different  sites  (Bizzi  et  al.  1991;  Giszter  et  al.  1993).  Next,  following  the 
experimental  procedure  described  in  (Giszter  et  al.  1993),  we  calculated  the  "active"  force  fields 
by  subtraction  the  resting  field  from  the  total  fields.  We  found  that  the  structure  of  force  fields  is 
similar  to  that  measured  experimentally  (Bizzi  et  al.  1991;  Giszter  et  al.  1993;  Mussa-Ivaldi  et  al. 
1994). 

Co-activation  of  distinct  sites.  Experimental  results  (Bizzi  et  al.  1991;  Mussa-Ivaldi  et  al. 
1994)  suggested  that  the  "vectorial  summation  rule"  is  valid  for  the  active  fields.  We  compared 
the  active  field  generated  by  the  simultaneous  activation  of  two  7N  units  to  the  vector  sum  of  the 
fields  generated  by  the  independent  activation  of  the  same  units.  As  it  was  suggested  by  Mussa- 
Ivaldi  et  al.  (1994),  to  quantify  the  similarity  between  the  sum  and  co-activation  fields,  we 
calculated  the  "cosine"  of  the  angle  between  these  two  sampled  fields,  and  used  for  this  measure 
the  value  of  0.90  as  a  threshold  for  similarity  (the  maximum  value  of  the  similarity  is  equal  to  1). 
There  are  four  units  at  the  IN  layer  and,  therefore,  there  are  six  different  pairs  of  them.  We  have 
found  that,  for  all  six  pairs,  the  similarity  between  the  sum  and  co-activation  fields  exceeded  the 
threshold  value  ranging  from  0.97  to  0.99. 

We  extended  this  analysis  comparing  the  vector  sum  of  fields  generated  by  two  random 
patterns  of  the  IN  activity  (the  activity  of  each  unit  is  a  continuous  variable  ranging  from  0  to  1) 
with  the  co-activation  field  generated  by  superposition  of  these  patterns.  We  calculated  the  sum 
and  co-activation  fields  for  10,000  pairs  of  random  patterns.  The  similarity  between  these  fields 
was  dependent  on  the  sampled  patterns  ranging  from  0.71  to  0.99.  Therefore,  the  co-activation 
field  generated  by  superposition  of  two  arbitrary  patterns  may  differ  dramatically  from  the  sum 
of  the  fields  generated  by  each  of  these  patterns  independently  (e.g.,  the  value  of  0.71 
corresponds  to  the  average  angle  between  sample  vectors  equal  to  45  degrees).  This  result  was 
expected  because  of  the  nonlinearity  of  the  model.  However,  the  number  of  cases  in  which  the 
similarity  was  below  the  threshold  value  0.90  comprised  less  than  15%  of  the  total  number  of 
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trials.  Averaged  over  all  trials  the  similarity  between  the  sum  and  co-activation  fields  was  0.96  + 
0.04,  which  is  consistent  with  the  average  similarity  across  the  entire  set  of  experimental  data 
0.94  +  0.05  (Mussa-Ivaldi  et  al.  1994). 

In  summary,  the  force  fields  have  been  calculated  within  the  same  model  that  has  been 
adjusted  and  tested  above.  None  of  the  parameters  has  been  changed,  neither  introduced 
additionally.  Nevertheless,  these  results  indicate  that  basic  features  of  the  force  fields  obtained  by 
the  spinal  stimulation  of  the  fi-og's  limb  (Bizzi  et  al.  1991;  Giszter  et  al.  1993;  Mussa-Ivaldi  et  al. 
1994),  such  as  a  convergence  of  total  fields  to  equilibrium  points,  typical  nonlinear  dependence 
of  the  force  direction  and  magnitude  upon  the  position  in  the  workspace,  and  the  "vectorial 
summation  aile",  are  well  reproduced  within  the  model. 

Conclusions 

The  complexity  of  real  biological  systems  controlling  motor  behaviors  is  such  that  no  practical 
models  attempt  to  describe  them  in  full  detail.  In  this  paper  we  focused  on  a  particular  aspect  of 
the  problem:  what  are  possible  mechanisms  underlying  the  integration  of  neuronal  commands 
firom  different  systems  to  produce  a  required  force.  We  aimed  to  develop  a  model  of  the 
"minimal  complexity"  that  would  include  those  factors  that  are  of  primarily  importance  for  the 
problem,  making  as  few  assumptions  as  possible  concerning  the  details  left  out  of  the  model. 

Although  the  degree  of  isomorphism  between  the  model  and  motor  physiology  should  not 
be  overstated  (for  example,  the  motoneuron  activation  function  (4)  does  not  contain  any  feedback 
component),  we  believe  that  the  model  described  in  this  paper  is  a  reasonable  compromise 
between  tractability  and  realism.  Indeed,  the  structure  of  the  model  is  really  simple:  it  consists  of 
a  standard  tliree-layered  feedforward  neural  network  connected  to  a  standard  two-joint-six- 
muscle  model  of  the  arm.  On  the  other  hand,  the  model  incorporates  the  following  biologically 
important  features:  force  signals  are  encoded  in  the  activity  of  directionally  tuned  neurons;  only  a 
few  distinct  sites  represent  the  level  of  spinal  cord  intemeurons;  synaptic  cormections  between 
neural  layers  result  in  converging  and  diverging  patterns  of  influence  (in  particular,  neuronal 
units  from  the  layer  representing  spinal  cord  intemeurons  make  synaptic  connections  with 
different  pools  of  motoneurons  and  activate  groups  of  muscles  in  a  weighted  fashion);  neuronal 
units  are  modelled  as  elements  with  nonlinear  activation  function;  the  construction  of  the  arm 
model  provides  a  nonlinear  dependence  of  elastic  properties  on  configuration  of  the  arm,  and 
muscles  are  modelled  by  nonlinear  springs. 

To  adjust  the  performance  of  the  model  we  used  two  sets  of  experimental  data.  One  set 
(the  stiffness  characteristics  of  human  arm)  was  used  to  train  the  network,  and  another  set  (the 
correlation  between  motor  cortical  activity  and  the  force  exerted  by  monkey,  and  the  stmcture  of 
restoring  force  fields  generated  by  stimulation  of  the  firog's  spinal  cord)  was  used  to  test  the 
model  performance.  We  have  shown  that  the  performance  of  the  model  with  a  fixed  set  of 
parameters  is  quantitatively  consistent  with  experimental  data  obtained  on  three  different 
organisms  using  different  experimental  paradigms  thus  suggesting  a  universal  principle  of 
organization  of  motor  control. 

The  model  predicts  a  specific  pattern  of  cormections  between  supraspinal  populations 
coding  forces  and  spinal  cord  intemeurons.  We  have  shown  that  a  supraspinal  signal  can  be 
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properly  translated  into  a  required  force  if  strengths  of  cortical-spinal  connections  are  weighted 
in  accordance  with  directional  preference  of  connected  units  (see  also  Georgopoulos,  1988  and 
Georgopoulos,  1994,  where  this  type  of  interactions  between  cortical  and  spinal  systems  was 
hypothesized).  In  this  paper,  the  cortical-spinal  connectivity  matrix  is  given  by  an  explicit 
relation  (10). 

In  fact,  our  investigation  originated  from  experimental  findings  (Georgopoulos  et  al. 

1992;  Mussa-Ivaldi  et  al.,  1994),  which  raised  a  question  of  how  essentially  nonlinear  system 
could  provide  a  linear  summation  of  forces  coding  by  different  neuronal  signals.  Obviously,  a 
nonlinear  system  cannot  carry  out  the  linear  summation  exactly,  merely  by  definition  of  what  the 
nonlinearity  is.  Indeed,  our  model,  which  incorporates  nonlinearities  at  all  levels  (1-11),  never 
performed  the  exact  vectorial  summation  of  forces.  However,  we  have  shown  that  adjustable 
parameters  of  the  model  (synaptic  connections)  can  be  chosen  in  a  way  that  the  force  generated 
by  the  sum  of  signals  is  nearly  equal  to  the  vector  sum  of  forces  generated  by  each  signal 
independently.  From  this  we  conclude  that  approximately  linear  smnmation  of  forces  observed  in 
the  experiments  (Georgopoulos  et  al.  1992;  Mussa-Ivaldi  et  al.  1994)  cannot  be  a  stringent  rule, 
but  may  reflect  a  prevalent  tendency,  which  could  be  provided  by  patterns  of  cortical-spinal  and 
spinal  intemeuron-motoneuronal  connections. 


2.B.  A  model  of  the  interactions  between  motor  cortical  and  spinal  networks  ' 

Abstract  ’ 

One  problem  in  motor  control  concerns  the  mechanism  whereby  the  central  nervous  system  ^ 
translates  the  motor  cortical  command  encoded  in  cell  activity  into  a  coordinated  contraction  of 
limb  muscles  to  generate  a  desired  motor  output.  This  problem  is  closely  related  to  the  design  of 
adaptive  systems  that  transform  neuronal  signals  chronically  recorded  from  the  motor  cortex  into 
the  physiologically  appropriate  motor  output  of  multijoint  prosthetic  limbs.  In  this  report,  we 
demonstrate  how  this  transformation  can  be  carried  out  by  an  artificial  neural  network  using  as 
command  signals  the  actual  impulse  activity  obtained  from  recordings  in  the  motor  cortex  of 
monkeys  during  the  performance  of  a  task  that  required  the  exertion  of  force  in  different 
directions.  The  network  receives  experimentally  measured  brain  signals  and  re-codes  them  into 
motor  actions  of  a  simulated  actuator  that  mimics  the  primate  arm.  The  actuator  responds  to  the 
motor  cortical  commands  with  surprising  fidelity  generating  forces  in  close  quantitative 
agreement  with  those  exerted  by  trained  monkeys,  in  both  the  temporal  and  spatial  domains. 
Moreover,  we  show  that  the  time-varying  motor  output  may  be  controlled  by  the  impulse  activity 
of  as  few  as  15  motor  cortical  cells.  These  results  outline  a  potentially  implementable 
computation  scheme  that  utilizes  raw  neuronal  signals  to  drive  artificial  mechanical  systems. 

Introduction 

Neurophysiological  studies  have  shown  that  motor  output  (movement  or  force)  produced  by  a 
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limb  is  related  to  the  motor  cortical  activity  and  that  the  population  activity  of  directionally  tuned 
cells  within  the  motor  cortex  can  be  used  to  predict  the  upcoming  motor  output  under  a  variety  of 
conditions.  However,  these  experimental  results  suggest  only  the  cortical  representation  of  the 
motor  command,  and  do  not  solve  the  problem  of  how  this  representation  could  be  implemented 
during  motor  output.  Here,  we  directly  address  the  question  of  how  the  cortical  neuronal  .signals 
could  be  used  to  drive  an  artificial  actuator  (e.g.,  a  prosthetic  limb)  such  that  the  motor  output  of 
actuator  would  correspond  to  the  performance  of  real  limb  driven  by  the  central  nervous  system. 
This  problem  is  made  difficult  by  the  fact  that  physiologically  usefiil  motor  output  requires  the 
coordinated  recruitment  of  multiple  degrees  of  freedom  distributed  over  artificial  limbs.  Below 
we  demonstrate  how  the  transformation  problem  can  be  solved  for  a  particular  example  of  motor 
action;  the  exertion  of  static  isometric  force. 

Methods 

Neuronal  signals:  The  motor  cortical  activity  used  in  this  study  as  command  signals 
came  from  single-cell  recordings  performed  previously  in  our  laboratory  for  other  purposes.  In 
that  experiment,  monkeys  were  trained  to  exert  2D  forces  on  an  isometric  handle  so  that  the  force 
developed  over  time  had  to  remain  in  the  visually  specified  direction  (the  instructed  direction) 
and  to  increase  in  magnitude  in  order  to  exceed  a  required  threshold.  Eight  instructed  directions 
uniformly  distributed  in  2D  space  were  employed.  In  the  present  study,  we  use  two  data  files. 
The  first  file  consisted  of  impulse  activity  of  75  directionally  tuned  cells  recorded  in  the  arm  area 
of  the  motor  cortex.  There  were  24  spike  trains  for  each  cell  (8  instructed  directions  of  3  trials 
each).  The  second  file  consisted  of  the  x-y  components  of  the  developed  force  measured  each  10 
ms  for  all  1800  trials  contained  in  the  first  file.  Therefore,  we  have  experimentally  measured 
motor  cortical  commands  and  resulting  motor  actions.  Our  aim  was  to  construct  a  model  that 
would  generate  the  same  motor  outputs  in  response  to  the  same  neuronal  commands. 

Transformation  algorithm:  We  exploit  a  processing  algorithm  developed  in  a  previous 
study.  The  impulse  activity  of  15  different  cells  is  presented.  These  spike  trains  were  recorded  in 
different  trials  but  for  the  same  instructed  direction  of  force  (-135°  in  2D  space).  The  cortical 
signal  must  be  transformed  into  the  motor  output  of  a  simulated  actuator.  For  example,  in  order 
to  be  consistent  with  the  experimental  data,  the  response  of  the  actuator  to  the  neuronal  command 
should  be  an  end-point  force  ‘exerted’  in  the  direction  -135°.  The  transformation  is  done  by  an 
artificial  neural  network  connected  to  the  actuator.  The  network  receives  experimentally 
measured  impulse  activity  as  a  time- varying  input  to  the  upper  layer  and  transforms  it  into  a 
time-varying  pattern  of  activity  at  the  output  layer  as  follows.  Each  neuron  at  the  intermediate 
layer  is  connected  to  all  neurons  at  the  upper  layer.  Intermediate  neurons  receive  the  input 
signals  as  the  sum  of  spikes  from  all  upper  neurons  weighted  with  corresponding  strengths  of 
synaptic  connections.  These  inputs  are  accumulated  over  the  time  and  transformed  into  activities 
of  intermediate  neurons  using  the  sigmoid  activation  function.  At  each  instant  of  time,  the 
collective  activity  of  intermediate  layer  produces  a  pattern  of  activity  at  the  output  layer  in 
accordance  with  the  connectivity  between  these  two  layers.  The  output  activity  of  the  network  is 
then  transformed  into  analog  commands.  Due  to  connections  between  the  output  units  and  the 
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actuator,  a  pattern  of  output  activity  generates  contractions  of  actuator  modeled  ‘muscles’  by 
means  of  changing  the  muscle  rest  lengths.  Let  the  actuator  be  at  the  equilibrium  state 
(equilibrium  joint  angles)  determined  by  an  initial  set  of  muscle  rest  lengths.  Now  suppose  that 
because  of  changes  in  the  network  output  activity  the  muscle  rest  lengths  have  changed.  If  the 
actuator  is  free  to  move,  then  a  new  equilibrium  state  would  eventually  be  reached.  Otherwise  (a 
possibility,  which  we  explore  here),  the  difference  between  the  previous  and  new  equilibrimn 
states  results  in  an  end-point  force,  which  we  refer  to  as  the  force  exerted  by  the  actuator  against 
an  immovable  object  (thick  arrow  at  the  end-point).  Given  a  set  of  new  rest  lengths,  the  direction 
and  magnitude  of  the  end-point  force  can  be  calculated  using  straightforward,  though 
cumbersome,  algebra.  In  the  framework  of  this  computational  scheme,  the  performance  of  the 
model  (i.e.,  the  relation  between  the  input  cortical  signals  and  the  forces  exerted)  depends  mainly 
on  the  network  connectivity,  which  must  provide  a  synergistic  activation  of  all  muscles  to 
generate  the  required  motor  output.  We  searched  the  corresponding  set  of  synaptic  weights  and 
activation  thresholds  using  the  algorithm  described  in  Lukashin  et  al.  (1996). 

Results 

After  the  connectivity  of  the  network  was  found  and  fixed,  we  tested  the  performance  of  the 
model  across  the  whole  set  of  experimental  data.  An  important  question  is  the  robustness  of  the 
computation  scheme  with  respect  to  the  variation  in  the  number  of  cells  generating  neuronal 
signals  {N)  and  with  respect  to  changes  in  the  cell  activity  from  trial  to  trial.  To  assess 
robustness,  we  adopted  an  approach  that  utilized  a  “bootstrapping”  procedme;  (i)  the  cells 
comprising  the  population  of  a  given  size  Wwere  selected  from  the  set  of  75  cells  (since  the 
imiformness  of  cell  preferred  directions  throughout  the  space  is  of  particular  importance  in 
reconstructing  the  directional  signal  encoded  in  the  population  activity,  all  cells  from  the  original 
data  set  were  divided  into  N  groups  in  accordance  with  their  preferred  directions  so  that  the  i-th 
group  would  contain  the  cells  with  preferred  direction  angles  ranging  from  27i(i-l)/A  to  2nilN\ 
to  comprise  the  population,  one  cell  from  each  group  was  drawn  randomly);  (ii)  for  each  drawn 
cell,  one  of  the  three  trials,  for  which  this  cell  was  recorded  from,  was  chosen  randomly;  (iii)  N 
spike  trains  corresponding  to  the  chosen  cells  and  trial  numbers  were  taken  from  the 
experimental  data  file  and  used  as  the  input  activity  of  the  network.  This  procedure  was  repeated 
1000  times  for  each  of  the  8  instructed  directions  using  new  sets  of  selected  cells  and  trial 
numbers.  The  size  of  the  population  A  varied  from  4  to  40.  To  evaluate  the  performance  of  the 
actuator  we  applied  the  same  criterion  as  in  the  experiment  (Georgopoulos  et  al.  1992):  the 
direction  of  force  developed  had  to  be  stable  and  not  differ  more  than  ±22.5*’  from  the  instructed 
direction.  We  have  found  (data  not  shown)  that  the  performance  improves  gradually  as  the  size 
of  population  A  increases,  revealing  a  tendency  for  saturation  as  A  approaches  15-20.  There  was 
a  close  similarity  between  the  actuator  performance  and  the  motor  action  of  the  real  arm.  We 
used  bootstraps  of  1000  randomly  selected  sets  of  cells  and  found  that  even  for  this  relatively 
small  number  of  command  cells  (A=  15)  the  mean  directions  of  forces  produced  by  the  actuator 
were  very  close  to  the  instructed  directions  (the  circular  correlation  coefficient  pj.  =  0.988,  n  =  8 
directions).  This  indicates  a  high  degree  robustness  of  the  transformation  algorithm. 
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Discussion 

In  a  sense,  the  control  of  isometric  force  chosen  in  the  present  study  represents  a  relatively 
simple  case,  if  compared  with  the  control  of  fast  movements.  For  the  latter  tasks,  the  problem  can 
be  additionally  complicated  by  the  necessity  of  compensating  for  the  dynamic  forces  (inertia, 
Coriolis  forces,  and  so  on)  acting  on  the  limb.  In  our  case,  there  is  no  motion  and,  therefore, 
there  is  no  need  to  allow  for  the  dynamic  forces.  Here,  we  use  a  computational  model  of  motor 
control  that  has  been  developed  previously  (Lukashin  et  al.  1996)  and  that  brings  together  a 
number  of  experimental  and  theoretical  observations  on  neuronal  activities  in  the  motor  cortex 
(Georgopoulos  et  al.  1993)  and  on  the  mechanics  of  multijoint  limb  control  (Mussa-Ivaldi  et  al. 
1985,  Bizzi  et  al.  1991,  Tsuji  et  al.  1995,  Bizzi  et  al.  1995).  In  our  previous  work  (Lukashin  et 
al.  1996)  we  concentrated  on  the  mechanisms  underlying  the  integration  of  neuronal  commands 
from  different  sources,  and  the  cortical  activity  was  modeled  by  artificial  idealized  signals.  The 
main  point  of  the  present  report  is  that  we  do  not  employ  the  artificial  signals.  Instead,  we  deal 
with  raw  experimental  data*  without  any  filtering,  approximation  or  pre-processing. 

Our  computation  scheme  successfully  decodes  these  signals,  and  the  time-vaiying  forces 
developed  by  the  actuator  are  very  similar  to  those  developed  by  the  monkey  arm.  Following  an 
initial  period  of  time,  which  lasts  100-200  ms,  the  direction  of  force  exerted  by  the  actuator 
stabilizes,  and  the  magnitude  of  force  increases.  The  stabilized  direction  of  force  is  close  to  the 
instructed  direction,  for  which  the  cortical  activity  was  recorded.  This  was  also  observed  for  the 
remaining  4  instructed  directions  for  this  particular  ensemble  of  cells  and  for  other  ensembles  of 
cells  {N>\S)  and  trial  numbers  that  were  selected  by  the  bootstrapping  procedure  as  described  in 
the  previous  section.  It  should  be  noted  that  different  patterns  of  cortical  activity  were  processed 
by  the  network  with  a  fixed  set  of  adaptive  parameters  (synaptic  weights  and  activation 
thresholds). 

To  evaluate  the  performance  of  the  model,  we  used  neuronal  activity  collected  from 
different  repeated  trials,  and  most  of  the  cells  were  not  recorded  simultaneously.  It  means  that 
spike  trains  of  different  cells  were  not  correlated  with  each  other.  For  real  time  applications  of 
the  transformation  scheme  suggested  here,  the  impulse  activity  of  different  cells  Avill  be  recorded 
simultaneously  during  an  individual  trial.  In  this  case,  different  cells  might  be  correlated,  and  the 
correlation  could  affect  the  performance  of  our  model.  Although  we  believe  that  the  correlation 
between  spike  trains  would  improve  the  performance,  this  is  still  a  conjecture  to  be  tested  in  the 
future.  Another  problem  that  may  arise  for  real  time  applications  is  that  some  of  motor  cortical 
cells  can  be  active  in  a  manner  consistent  with  the  incoming  motor  action  during,  for  example, 
instructed  delays.  In  such  a  situation,  the  motor  cortical  activity  alone  might  not  be  sufficient  to 
generate  the  desired  time-course  of  motor  actions,  and  an  additional  signal  that  triggers  the 
system  is  needed. 

Conclusion 

If  it  becomes  practical  to  chronically  record  motor  cortical  activity  as  command  signals 
for  electronically  driven  prostheses  or  other  mechanical  systems,  then  the  transformation  of  such 
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signals  into  appropriate  motor  output  would  become  an  important  issue.  The  reliable  decoding 
and  transformation  of  neuronal  ensemble  activity  recorded  in  behaving  animals  as  demonstrated 
here  suggests  that  the  use  of  biologically  inspired  neural  networks  to  transform  raw  cortical 
signals  into  the  motor  output  of  a  multijoint  artificial  limb  is  both  feasible  and  practical. 
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